# Copyright (C) 2022 VRVis.
# All rights reserved.
# Contact: VRVis Forschungs-GmbH (office@vrvis.at)
# 
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# 1. Redistributions of source code must retain the above copyright
#    notice, this list of conditions and the following disclaimer.
# 2. Redistributions in binary form must reproduce the above copyright
#    notice, this list of conditions and the following disclaimer in the
#    documentation and/or other materials provided with the distribution.
# 3. All advertising materials mentioning features or use of this software
#    must display the following acknowledgement:
#    This product includes software developed by the VRVis Forschungs-GmbH.
# 4. Neither the name of the VRVis Forschungs-GmbH nor the
#    names of its contributors may be used to endorse or promote products
#    derived from this software without specific prior written permission.
#
#setup

workingDir<-"C:/Users/ganglberger/Documents/Evolution/evo_paper_v3_olga_data/R_MAIN_CODE"
knitr::opts_knit$set(root.dir = normalizePath(workingDir)) 

setwd(workingDir)

set.seed(42)
#load libraries

#Runs under R 4.0.3
options(java.parameters = "- Xmx2024m")
library("xlsx")
library("corrplot")
## corrplot 0.84 loaded
library("R.matlab")
## R.matlab v3.6.2 (2018-09-26) successfully loaded. See ?R.matlab for help.
## 
## Attaching package: 'R.matlab'
## The following objects are masked from 'package:base':
## 
##     getOption, isOpen
library("rjson")
library("gplots")
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
library("png")
library("e1071")
library("rsvg")
library("XML")
library("igraph")
## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
library('knitr')
library("biomaRt")
library('org.Hs.eg.db')
## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:igraph':
## 
##     normalize, path, union
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which, which.max, which.min
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: IRanges
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:gplots':
## 
##     space
## The following object is masked from 'package:base':
## 
##     expand.grid
## 
## Attaching package: 'IRanges'
## The following object is masked from 'package:grDevices':
## 
##     windows
## 
#Read DNDS-ratio data

dnds<-read.csv2(dataFile,stringsAsFactors=FALSE,dec = ".",sep=",",header=length(colsToUseNames)==0)
# values below zero mean:
#-1 - unstable dn/ds in 5 runs
#-2 - short edge, unstable
#-3 - 0 values, dN=0
#-4 - 999 values, dS=0
#-5 identical sequences (the majority between HUMAN and NEANDERTHAL)

if(length(colsToUseNames)>0){
  colnames(dnds)[colsToUse]<-colsToUseNames
}

humanEntrez<-dnds[,2]

print(paste0("Dimension of data: ",dim(dnds)[1]," x ",dim(dnds)[2]))
## [1] "Dimension of data: 8978 x 25"
head(dnds)
#Read ontology (hierarchy of regions + a color palette for the output of MRI style plots)

documentHuman <- fromJSON(file="../storage/ontologyHuman.json", method='C')
#gets acronym of a region given by a region id (from ontology file)
getAcronymByID <- function(child,ID){
  if(!is.null(child)){
    if(child$id==ID){
      return(child)
    }
  }
  
  for(actchild in child$children){
    res<-getAcronymByID(actchild,ID)
    if(!is.null(res)){
      if(res$id==ID){
        return(res)
      }
    }
  }
  return(NULL)
}

getIDbyName <- function(child,name){
  if(!is.null(child)){
    if(child$name==name){
      return(child)
    }
  }
  
  for(actchild in child$children){
    res<-getIDbyName(actchild,name)
    if(!is.null(res)){
      if(res$name==name){
        return(res)
      }
    }
  }
  return(NULL)
}

jet.colors <-
  colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                     "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
#Read Human data and filter it for genes that are also in the DNDS-ratio file

#gives a vector with the size of the atlas (1D vector of human brain atlas biopsy sites) where every position
#is zero, execpt for positions where the region (given by the ID) can be found
getAtlasRegionsOfIDHuman <- function(atlasRegions,ID){
  newAtlasRegions<-atlasRegions==ID
  
  childrens<-getAcronymByID(documentHuman$msg[[1]],ID)$children
  while(length(childrens)>0){
    newChildren<-c()
    for(i in 1:length(childrens)){
      newAtlasRegions<-newAtlasRegions|(atlasRegions==childrens[[i]]$id)
      newChildren<-c(newChildren,getAcronymByID(documentHuman$msg[[1]],childrens[[i]]$id)$children)
    }
    
    childrens<-newChildren
  }
  return(newAtlasRegions)
}

atlasRegionsHuman<-readMat("../storage/atlasRegionsHuman.mat") #1D vector of human brain atlas biopsy sites
regionInfo<-cbind(atlasRegionsHuman$indexX,atlasRegionsHuman$indexY,atlasRegionsHuman$indexZ)
atlasRegionsHuman<-as.vector(atlasRegionsHuman$atlasRegions)

regionName<-c()
regionShortName<-c()

for(rid in atlasRegionsHuman){
  a<-getAcronymByID(documentHuman$msg[[1]],rid)
  regionName<-c(regionName,a$name)
  regionShortName<-c(regionShortName,a$acronym)
}

regionInfo<-cbind(regionShortName,regionName,regionInfo)
colnames(regionInfo)<-c("acronym","name","mniX","mniY","mniZ")



if(!exists("expressionHuman") && !exists("entrezIDsOfRowsHuman")){ #this if just makes it faster to reload biopsy site level gene expression data, since readMat is quite slow
  all_genes<-readMat("../storage/all_genes_expression_Human.mat")
  expressionHuman<-all_genes$expressionMatrix
  entrezIDsOfRowsHuman<-all_genes$entrezIDsOfRows
  
  remove(all_genes)
  gc()
  
  
  expressionHuman<-t(apply(expressionHuman,1,function(x){(x-median(x,na.rm=TRUE))/mad(x,na.rm=TRUE)})) #normalize robust per biopsy site
  expressionHuman<-apply(expressionHuman,2,function(x){(x-median(x,na.rm=TRUE))/mad(x,na.rm=TRUE)}) #normalize robust per gene (i.e. over the brain)
}
## 
## Attaching package: 'Matrix'
## The following object is masked from 'package:S4Vectors':
## 
##     expand
print(paste0("Dimension of human data: ",dim(expressionHuman)[1]," x ",dim(expressionHuman)[2]))
## [1] "Dimension of human data: 3702 x 20787"
entrezIDsInAMBAandDNDSHuman<-intersect(unique(humanEntrez),unique(entrezIDsOfRowsHuman))
print(paste0("Human genes that are in ABA and DNDS-ratio file: ",length(entrezIDsInAMBAandDNDSHuman)))
## [1] "Human genes that are in ABA and DNDS-ratio file: 8977"
#Filter genes that are human (and are in the DNDS ratio file)

humanEntrezInAmba<-unique(humanEntrez[is.element(humanEntrez,entrezIDsInAMBAandDNDSHuman)])
print(paste0("Genes that are in human (and are in the DNDS ratio file): ",length(humanEntrezInAmba), " unique: ",length(unique(humanEntrezInAmba))))
## [1] "Genes that are in human (and are in the DNDS ratio file): 8977 unique: 8977"
#Select columns (branches) that should be used. Take only rows (genes) that have no NAs in the branches and are not all 0
#Take also the mean for genes that are not unique for a human

dnds[,colsToUse]<-apply(dnds[,colsToUse],2,function(x){as.numeric(x)})
dnds_org<-dnds

dnds[,colsToUse][(dnds[,colsToUse]==(-3)) & !is.na(dnds[,colsToUse])]<-0 #set DNDS ratios wih a -3 value (dN=0, dS>0) to 0

dnds[,colsToUse][dnds[,colsToUse]<0 & !is.na(dnds[,colsToUse])]<-NA # set other special values below 0 to NA
dnds_org[,colsToUse][dnds_org[,colsToUse]<0 & !is.na(dnds_org[,colsToUse])]<-NA


print(paste0("Use cols: "))
## [1] "Use cols: "
print(paste0(names(dnds)[colsToUse]))
##  [1] "MOUSE"       "OTOGA"       "CALJA"       "MACMU"       "NOMLE"      
##  [6] "GORGO"       "PANTR"       "DENISOVAN"   "NEANDERTHAL" "HUMAN"      
## [11] "MO_CMNGPDNH" "MOC_MNGPDNH" "MOCM_NGPDNH" "MOCMN_GPDNH" "MOCMNG_HNDP"
## [16] "MOCMNGP_HND" "MOCMNGPD_HN" "MOCMNP_HNDG" "MOCMNDHN_PG" "MOCMNGPH_ND"
## [21] "MOCMNGPN_HD"
print(paste0("Rows before filter: ",dim(dnds)[1]))
## [1] "Rows before filter: 8978"
dnds_per_gene<-matrix(0,nrow=length(humanEntrezInAmba),ncol=length(colsToUse))
for(actCol in 1:length(colsToUse)){
  for(actRow in 1:length(humanEntrezInAmba)){
    dnds_per_gene[actRow,actCol]<-mean(dnds[humanEntrez==humanEntrezInAmba[actRow],colsToUse[actCol]],na.rm=TRUE)
  }
}

dnds_per_gene[is.nan(dnds_per_gene)]<-NA

print(paste0("Rows after filter for unique human of the selected columns: ",dim(dnds_per_gene)[1]))
## [1] "Rows after filter for unique human of the selected columns: 8977"
rownames(dnds_per_gene)<-humanEntrezInAmba
colnames(dnds_per_gene)<-names(dnds)[colsToUse]


print(paste0("Dimension of data: ",dim(dnds_per_gene)[1]," x ",dim(dnds_per_gene)[2]," length of entrez (should be first dimension: ",length(humanEntrezInAmba)))
## [1] "Dimension of data: 8977 x 21 length of entrez (should be first dimension: 8977"
print("Amount of rows not NA")
## [1] "Amount of rows not NA"
apply(as.matrix(dnds_per_gene),2,function(x){sum(!is.na(x))})
##       MOUSE       OTOGA       CALJA       MACMU       NOMLE       GORGO 
##        8936        8917        8865        8832        8737        7875 
##       PANTR   DENISOVAN NEANDERTHAL       HUMAN MO_CMNGPDNH MOC_MNGPDNH 
##        7261        8836        1097        2222        8786        8350 
## MOCM_NGPDNH MOCMN_GPDNH MOCMNG_HNDP MOCMNGP_HND MOCMNGPD_HN MOCMNP_HNDG 
##        8013        8149        3899        6893        1199         639 
## MOCMNDHN_PG MOCMNGPH_ND MOCMNGPN_HD 
##         652        1507          54
#Rank-Normalize dnds column wise

print(paste0("Column means/sd before normalization: "))
## [1] "Column means/sd before normalization: "
cbind(apply(as.matrix(dnds_per_gene),2,function(x){mean(x,na.rm=TRUE)}),apply(as.matrix(dnds_per_gene),2,function(x){sd(x,na.rm=TRUE)}))
##                  [,1]        [,2]
## MOUSE       0.1253018  0.11636026
## OTOGA       0.1619856  0.16315102
## CALJA       0.2087074  0.23122606
## MACMU       0.1904524  0.24372696
## NOMLE       0.2347050  0.32149102
## GORGO       0.2540529  0.36566381
## PANTR       0.2336738  0.34285142
## DENISOVAN   0.5903055  0.48754410
## NEANDERTHAL 0.1045043  0.23628585
## HUMAN       0.1066694  0.22554817
## MO_CMNGPDNH 0.1745729  0.21954930
## MOC_MNGPDNH 0.2358727  0.63212933
## MOCM_NGPDNH 0.2188350  0.48567944
## MOCMN_GPDNH 0.2133775  0.40029925
## MOCMNG_HNDP 0.4113882 15.96900842
## MOCMNGP_HND 0.2289946  0.35842950
## MOCMNGPD_HN 0.0923868  0.22637260
## MOCMNP_HNDG 0.1120567  0.62693088
## MOCMNDHN_PG 0.0984602  0.23519968
## MOCMNGPH_ND 0.1018453  0.21166113
## MOCMNGPN_HD 0.0161238  0.06998989
normalized_dnds_per_gene<-(apply(as.matrix(dnds_per_gene),2,function(x){
  y<-x[!is.na(x)]
  y[y>0]<-rank(y[y>0],ties.method="min",na.last="keep")/sum(!is.na(y[y>0]))
  x[!is.na(x)]<-y
  x
  }))
normalized_dnds_per_gene[dnds_per_gene==0]<-0 
colnames(normalized_dnds_per_gene)<-colnames(dnds_per_gene)



print(paste0("Column means/sd after normalization: "))
## [1] "Column means/sd after normalization: "
cbind(apply(as.matrix(normalized_dnds_per_gene),2,function(x){mean(x,na.rm=TRUE)}),apply(as.matrix(normalized_dnds_per_gene),2,function(x){sd(x,na.rm=TRUE)}),apply(as.matrix(normalized_dnds_per_gene),2,function(x){sum(x==min(x,na.rm=TRUE),na.rm=TRUE)}))
##                   [,1]      [,2] [,3]
## MOUSE       0.48897649 0.2948124  198
## OTOGA       0.48026180 0.2992517  353
## CALJA       0.46373308 0.3068162  644
## MACMU       0.42912081 0.3193213 1253
## NOMLE       0.41644684 0.3228531 1461
## GORGO       0.35695197 0.3325483 2254
## PANTR       0.31400591 0.3328298 2702
## DENISOVAN   0.49456703 0.2917929   97
## NEANDERTHAL 0.13445761 0.2679131  803
## HUMAN       0.13883888 0.2709081 1606
## MO_CMNGPDNH 0.45788714 0.3092364  741
## MOC_MNGPDNH 0.38359235 0.3295723 1945
## MOCM_NGPDNH 0.32852826 0.3333506 2749
## MOCMN_GPDNH 0.33770972 0.3333556 2646
## MOCMNG_HNDP 0.20376987 0.3072317 2311
## MOCMNGP_HND 0.29878067 0.3315979 2775
## MOCMNGPD_HN 0.11509591 0.2523414  924
## MOCMNP_HNDG 0.10563380 0.2441723  505
## MOCMNDHN_PG 0.12269939 0.2590456  493
## MOCMNGPH_ND 0.13171697 0.2657447 1111
## MOCMNGPN_HD 0.03703704 0.1671900   51
#Plot and write not normalized data

boxplot(dnds_per_gene,main="Not normalized")

boxplot(dnds_per_gene,outline=FALSE,main="Not normalized (without outlier)")

png(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_boxplot.png"),width=round(1200*0.45),height=300)
boxplot(dnds_per_gene,outline=FALSE,main="Not normalized (without outlier)")
dev.off();
## png 
##   2
write.csv2(dnds_per_gene,paste0('../storage/',paste0("csvs/",prefix),"dnds_per_gene.csv"))
#Plot and write not normalized data with quantiles

par(mar=c(8,5,1,1))
boxplot(dnds_per_gene,main="Not normalized",las=2)

boxplot(dnds_per_gene,outline=FALSE,main="Not normalized (without outlier)",las=2)

n <- ncol(dnds_per_gene)
# width of each boxplot is 0.8
x0s <- 1:n - 0.4
x1s <- 1:n + 0.4
# these are the y-coordinates for the horizontal lines
# that you need to set to the desired values.
y0s <- apply(dnds_per_gene,2,function(x){quantile(x,0.9,na.rm=TRUE)})


# add segments
segments(x0 = x0s, x1 = x1s, y0 = y0s, col = "red")

png(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_boxplot.png"),width=round(1200*0.45),height=300)
par(mar=c(8,5,1,1))
boxplot(dnds_per_gene,outline=FALSE,main="Not normalized (without outlier)",las=2)
dev.off();
## png 
##   2
png(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_boxplot_90percentile.png"),width=round(1200*0.45),height=300)
par(mar=c(8,5,1,1))
boxplot(dnds_per_gene,outline=FALSE,main="Not normalized (without outlier)",las=2)
segments(x0 = x0s, x1 = x1s, y0 = y0s, col = "red")
dev.off();
## png 
##   2
#Plot and write normalized data

boxplot(normalized_dnds_per_gene,outline=FALSE,main="normalized (without outlier)")

boxplot(normalized_dnds_per_gene,outline=TRUE,main="normalized (with outlier)")

write.csv2(normalized_dnds_per_gene,paste0('../storage/',paste0("csvs/",prefix),"normalized_dnds_per_gene.csv"))
#Plot DNDS gene wide heatmap

normalized_dnds_per_gene_with_reordered<-normalized_dnds_per_gene[,c(1:19,21,20)]

heatmap.2(normalized_dnds_per_gene_with_reordered,
          density.info="none",
          trace="none",
          scale="none",
          dendrogram="row",
          na.color="#F0F0F0",
          Colv=NA,
          #keysize = 0.8,
          #key.par = list(mar=c(4,0.3,8,0.3),cex=6),
          key.title=NA,
          key.xlab="Ranked DN/DS",
          cexCol=1,
          srtCol= 90,
          labRow = FALSE,
          col=c(colorpanel(100,"#ffffb2","#feb24c","#bd0026")),
          margins = c(20,5),
          distfun = function(mat) {
              edist <- dist(mat)
              edist[which(is.na(edist))] <- max(edist, na.rm=TRUE) * 1.1 
              return(edist)
          },
          adjCol = c(1,1))

#just also save it as png or svg
#png(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_heatmap.png"),width=round(3200*0.45),height=3200)
svg(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_heatmap.svg"),width=9,height=20)
heatmap.2(normalized_dnds_per_gene_with_reordered,
          density.info="none",
          trace="none",
          scale="none",
          dendrogram="row",
          na.color="#F0F0F0",
          Colv=NA,
          #keysize = 0.8,
          #key.par = list(mar=c(4,0.3,8,0.3),cex=6),
          key.title=NA,
          key.xlab="Ranked DN/DS",
          cexCol=1,
          srtCol= 90,
          labRow = FALSE,
          col=c(colorpanel(100,"#ffffb2","#feb24c","#bd0026")),
          margins = c(20,5),
          distfun = function(mat) {
              edist <- dist(mat)
              edist[which(is.na(edist))] <- max(edist, na.rm=TRUE) * 1.1 
              return(edist)
          },
          adjCol = c(1,1))
dev.off()
## png 
##   2
dnds_per_gene_with_reordered<-dnds_per_gene[,c(1:19,21,20)]
dnds_per_gene_with_reordered[dnds_per_gene_with_reordered>1.5]<-1.5
heatmap.2(dnds_per_gene_with_reordered,
          density.info="none",
          trace="none",
          scale="none",
          dendrogram="row",
          na.color="#F0F0F0",
          Colv=NA,
          #keysize = 0.8,
          #key.par = list(mar=c(4,0.3,8,0.3),cex=6),
          key.title=NA,
          key.xlab="Ranked DN/DS",
          cexCol=1,
          srtCol= 90,
          labRow = FALSE,
          col=c(colorpanel(100,"#ffffb2","#feb24c","#bd0026")),
          margins = c(20,5),
          distfun = function(mat) {
              edist <- dist(mat)
              edist[which(is.na(edist))] <- max(edist, na.rm=TRUE) * 1.1 
              return(edist)
          },
          adjCol = c(1,1))

#just also save it as png or svg
#png(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_heatmap_raw.png"),width=round(3200*0.45),height=3200)
svg(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_heatmap_raw.svg"),width=9,height=20)
heatmap.2(dnds_per_gene_with_reordered,
          density.info="none",
          trace="none",
          scale="none",
          dendrogram="row",
          na.color="#F0F0F0",
          Colv=NA,
          #keysize = 0.8,
          #key.par = list(mar=c(4,0.3,8,0.3),cex=6),
          key.title=NA,
          key.xlab="Ranked DN/DS",
          cexCol=1,
          srtCol= 90,
          labRow = FALSE,
          col=c(colorpanel(100,"#ffffb2","#feb24c","#bd0026")),
          margins = c(20,5),
          distfun = function(mat) {
              edist <- dist(mat)
              edist[which(is.na(edist))] <- max(edist, na.rm=TRUE) * 1.1 
              return(edist)
          },
          adjCol = c(1,1))
dev.off()
## png 
##   2
#get names of genes
humanGeneNames<-unlist(mapIds(org.Hs.eg.db, paste(humanEntrezInAmba), 'SYMBOL', 'ENTREZID'))
## 'select()' returned 1:1 mapping between keys and columns
write.csv2(humanEntrezInAmba,paste0('../storage/',paste0("csvs/",prefix),"humanEntrezIDs.csv"))
#Generate genes lists for predicting functional maps WEIGHTED

addEntrezIDsHumanWeighted <- function(humanEntrezInAmba,gene_list_name_entrez_human,ratios){
  setEntrezHuman<-c()

  for(entrezIndex in (1:length(humanEntrezInAmba))){
    humanEntrez<-humanEntrezInAmba[entrezIndex]
  
    if(!is.na(humanEntrez)){
      if(length(setEntrezHuman)<length(unique(c(setEntrezHuman,humanEntrez)))){
        setEntrezHuman<-c(setEntrezHuman,humanEntrez)
        geneName<-unique(humanGeneNames[humanEntrezInAmba==humanEntrez])
        geneName<-geneName[!is.na(geneName)]
        if(length(geneName)==0){
          geneName<-paste0("entrez_",humanEntrez)
        }else{
          geneName<-geneName[1]
        }
        
        #humanName<-strsplit(geneName,"_")[[1]][1]
        humanName<-geneName
        humanName<-gsub("\\-", "_",humanName)
        humanName<-gsub("\\.", "_",humanName)
        humanName<-gsub("\\ ", "_",humanName)
        if(!is.na(ratios[entrezIndex]) && (ratios[entrezIndex]>0)){ #exclude genes with 0 weightings (just makes the list smaller and therefore the algorithm more perfomant)
          gene_list_name_entrez_human<-rbind(gene_list_name_entrez_human,c(humanName,humanEntrez,ratios[entrezIndex]))
        }
      }
    }
    
  
  }
  return(gene_list_name_entrez_human)
}

gene_list_name_entrez_human<-matrix("",0,3)

set.seed(42)

for(branch in transistionsToUseForFunctionalMap){
  print(paste0("branch: ",branch))
  print(paste0("GeneAmount:",sum(!is.na(normalized_dnds_per_gene[,branch]) & (normalized_dnds_per_gene[,branch]>0))))
  gene_list_name_entrez_human<-rbind(gene_list_name_entrez_human,c(paste0("dndsweighted_trans",branch),"",""))

  gene_list_name_entrez_human<-addEntrezIDsHumanWeighted(humanEntrezInAmba,gene_list_name_entrez_human,(normalized_dnds_per_gene[,branch]))
gene_list_name_entrez_human<-rbind(gene_list_name_entrez_human,c("","",""))

}
## [1] "branch: 11"
## [1] "GeneAmount:8045"
## [1] "branch: 12"
## [1] "GeneAmount:6405"
## [1] "branch: 13"
## [1] "GeneAmount:5264"
## [1] "branch: 14"
## [1] "GeneAmount:5503"
## [1] "branch: 15"
## [1] "GeneAmount:1588"
## [1] "branch: 16"
## [1] "GeneAmount:4118"
## [1] "branch: 17"
## [1] "GeneAmount:275"
## [1] "branch: 10"
## [1] "GeneAmount:616"
dir.create("../storage/setsToRunHuman")
## Warning in dir.create("../storage/setsToRunHuman"): '..\storage\setsToRunHuman'
## already exists
write.table(gene_list_name_entrez_human,paste0(paste0("../storage/setsToRunHuman",'/',strsplit(functionalMapsHuman,"/")[[1]][length(strsplit(functionalMapsHuman,"/")[[1]])],".csv")),sep=";", row.names = FALSE, col.names = FALSE, quote = FALSE)
#Load task fmri data

taskfmri<-readMat(paste0("../storage/taskfMRI.mat"))$taskfMRI
tasknames<-unlist(readMat(paste0("../storage/taskfMRI.mat"))$tasknames)
tasknames <- gsub("tfMRI_", "", tasknames)

#Remove the contrast tasks i.e. only take the direct tasks
dif_tasks<-c(3,6,9,10,17:22,25,28,39,44:47)
tasknames<-tasknames[-dif_tasks]
taskfmri<-taskfmri[,-dif_tasks]
colnames(taskfmri)<-tasknames

#define taskgroups
colnames(taskfmri)[5:6] <- paste("LA-AVG", colnames(taskfmri)[5:6], sep = "_")
colnames(taskfmri)[c(7,9)] <- paste("MO-RF+LF", colnames(taskfmri)[c(7,9)], sep = "_")
colnames(taskfmri)[11] <- paste("MO-T", colnames(taskfmri)[11], sep = "_")
colnames(taskfmri)[c(8,10)] <- paste("MO-RH+LH", colnames(taskfmri)[c(8,10)], sep = "_")
colnames(taskfmri)[12] <- paste("MO-RH+LH", colnames(taskfmri)[12], sep = "_")
colnames(taskfmri)[16] <- paste("SO-TOM", colnames(taskfmri)[16], sep = "_")
colnames(taskfmri)[1] <- paste("EM-FACES", colnames(taskfmri)[1], sep = "_")
colnames(taskfmri)[2] <- paste("EM-SHAPE", colnames(taskfmri)[2], sep = "_")
taskfmri[,3]<-apply(taskfmri[,3:4],1,mean)
colnames(taskfmri)[4] <- paste("GA-REWARD", colnames(taskfmri)[4], sep = "_")
colnames(taskfmri)[3] <- paste("GA-AVG", "GAMBLING_AVG", sep = "_")
colnames(taskfmri)[13:14] <- paste("RE-AVG", colnames(taskfmri)[13:14], sep = "_")
colnames(taskfmri)[c(18,22,28)] <- paste("WM-FACE", colnames(taskfmri)[c(18,22,28)], sep = "_")
colnames(taskfmri)[c(17,19:21,23:27,29:30)] <- paste("WM-REST", colnames(taskfmri)[c(17,19:21,23:27,29:30)], sep = "_")

#Remove tasks that were not used in the paper
taskfmri<-taskfmri[,-15]
taskfmri<-taskfmri[,-2]

taskfmri[is.na(taskfmri)]<-0
taskfmri[abs(taskfmri)<5]<-0 #removing insiginifcant values (smaller than 5 sigmas)

heatmap.2(cor(taskfmri),cexCol=0.5,cexRow=0.5,scale="none",density.info="none", trace="none")

colorsOfBiopsyHuman<-c()
for(actRegion in 1:length(atlasRegionsHuman)){
  colorsOfBiopsyHuman<-c(colorsOfBiopsyHuman,paste0("#",getAcronymByID(documentHuman$msg[[1]],atlasRegionsHuman[[actRegion]][1])$color_hex_triplet))
}

namesOfBiopsyHuman<-c()
for(actRegion in 1:length(atlasRegionsHuman)){
  namesOfBiopsyHuman<-c(namesOfBiopsyHuman,paste0(getAcronymByID(documentHuman$msg[[1]],atlasRegionsHuman[[actRegion]][1])$name))
}

taskClasses<-sapply(colnames(taskfmri),function(x){sapply(strsplit(x, "_"),"[[",1)})

taskGroupfmri<-sapply(unique(taskClasses),function(x){
  if(sum(x==taskClasses)>1){
    rowMeans(taskfmri[,x==taskClasses])
    }else{
      taskfmri[,x==taskClasses]
    }})

rownames(taskGroupfmri)<-namesOfBiopsyHuman
 
# heatmap.2(t(taskGroupfmri[order(colorsOfBiopsyHuman),]),hclustfun=function(x){hclust(x,method="ward.D2")},main="taskGroupFMRI cortical regions",scale="none",density.info="none", trace="none",ColSideColors = sort(colorsOfBiopsyHuman),Colv=NA,col=colorpanel(100, "blue","white", "red"),margins=c(6,10),cexRow=0.8,cexCol=0.1,dendrogram="row")
#  
# heatmap.2(t(abs(taskGroupfmri)[order(colorsOfBiopsyHuman),]),hclustfun=function(x){hclust(x,method="ward.D2")},main="taskGroupFMRI cortical regions",scale="none",density.info="none", trace="none",ColSideColors = sort(colorsOfBiopsyHuman),Colv=NA,col=colorpanel(100, "white","orange", "red"),margins=c(6,10),cexRow=0.8,cexCol=0.1,dendrogram="row")
#Load network data

#Creating networks from literature
networks<-matrix(0,nrow=sum(atlasRegionsHuman>0),ncol=11)
colnames(networks)<-c("cortico-limbic_network","mPFC-AMY_network","mPFC-Acb_network","deafult_mode_network","salience_network","central_executive_network","sensorimotor_network","visuospatial_network","fronto-parietal","dorsal_attention","ventral_attention")

networksIDHumanList<-list()

networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4035,4021,4050,4028,4268,4220,4249,9391,4212,4205,4156,4096,4111,4118,4071,4125,4278,4287,4293,4392,4140,4147,4053,4327)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4021,4050,4047,4897,4062,4327)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4021,4050,4047,4897,4062,4290)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4021,4050,4047,4897,4062,4235,4228,4242,4147,4111,4249)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4225,4272,4290,9066,9072,4540,12927,4320,4313)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4021,4028,4085,9372,4278,12920)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4021,4272,4010,4085)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4028,4085,4181,4118,4184)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4028,4096,4103)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4104,4096,4028)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4035,4133,4103)


for(actNetwork in 1:length(networksIDHumanList)){
  networks[,actNetwork]<-rep(FALSE,length(atlasRegionsHuman))
  for(actSubRegionID in networksIDHumanList[[actNetwork]]){
    networks[,actNetwork]<-networks[,actNetwork]|(getAtlasRegionsOfIDHuman(atlasRegionsHuman,actSubRegionID)>0)
  }
}

#heatmap.2(t(abs(networks)[order(colorsOfBiopsyHuman),]),hclustfun=function(x){hclust(x,method="ward.D2")},main="literature network regions",scale="none",density.info="none", trace="none",ColSideColors = sort(colorsOfBiopsyHuman),Colv=NA,col=colorpanel(100, "white","orange", "red"),margins=c(6,10),cexRow=0.8,cexCol=0.1,dendrogram="row")
#Visualize branches to networks and tasks combined on region level

regionsIDsHumanLarge<-list()
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4035)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4021)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4050)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4028)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4268)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4220)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4249)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4212)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4205)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4156)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4096)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4111)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4118)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4071)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4125)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4140)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4147)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4053)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4047)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4897)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4062)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4235)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4228)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4225)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4085)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4010)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4181)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4184)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4278)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4287)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4293)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4392)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4327)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4290)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(9066)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(9072)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4540)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4313)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4321)

namesOfRegionsLarge<-c("IFG",
"SFG",
"MOrG",
"MFG",
"Ins",
"CgG",
"HiF",
"SOG",
"IOG",
"FuG",
"SPL",
"AnG",
"Pcu",
"PCLa",
"PCLp",
"MTG",
"ITG",
"AOrG",
"GRe",
"SRoG",
"SCG",
"CgGr",
"CgGp",
"CgGf",
"PoG",
"PrG",
"OP",
"Cun",
"Cd",
"Pu",
"GP",
"TH",
"Amg",
"Acb",
"VTA",
"SN",
"Hy",
"BSTl",
"Cl")

colorsOfRegionsLarge<-c()

combinedNetworks<-cbind(abs(taskGroupfmri),abs(networks))
combinedNetworksRegionWise<-matrix(0,0,ncol(combinedNetworks))

for(actRegionIndex in 1:length(regionsIDsHumanLarge)){
  actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[actRegionIndex][1])>0
  a<-getAcronymByID(documentHuman$msg[[1]],regionsIDsHumanLarge[actRegionIndex][1])
  colorsOfRegionsLarge<-c(colorsOfRegionsLarge,paste0("#",a$color_hex_triplet))
  combinedNetworksRegionWise<-rbind(combinedNetworksRegionWise,apply(combinedNetworks,2,function(x){
    max(x[actRegionAtlas])
  }))
}
rownames(combinedNetworksRegionWise)<-namesOfRegionsLarge


write.csv2(t(combinedNetworksRegionWise[order(colorsOfRegionsLarge),1:ncol(abs(taskGroupfmri))]),paste0('../storage/',paste0("csvs/"),"task_network_regions.csv"))
heatmap.2(t(combinedNetworksRegionWise[order(colorsOfRegionsLarge),1:ncol(abs(taskGroupfmri))]),hclustfun=function(x){hclust(x,method="ward.D2")},main="task network regions",scale="none",density.info="none", trace="none",ColSideColors = sort(colorsOfRegionsLarge),Colv=NA,col=colorpanel(100, "white","orange", "red"),margins=c(6,10),cexRow=0.8,cexCol=0.8,dendrogram="row")

write.csv2(t(combinedNetworksRegionWise[order(colorsOfRegionsLarge),(1:ncol(abs(taskGroupfmri)))+ncol(abs(networks))]),paste0('../storage/',paste0("csvs/"),"literuature_network_regions.csv"))
heatmap.2(t(combinedNetworksRegionWise[order(colorsOfRegionsLarge),(1:ncol(abs(taskGroupfmri)))+ncol(abs(networks))]),hclustfun=function(x){hclust(x,method="ward.D2")},main="literature network regions",scale="none",density.info="none", trace="none",ColSideColors = sort(colorsOfRegionsLarge),Colv=NA,col=colorpanel(100, "white","orange", "red"),margins=c(6,10),cexRow=0.8,cexCol=0.8,dendrogram="row")

#Create gene-wise correlation to tasks and networks for biclustering

  resultsTableCorrelation<-matrix(0,length(humanEntrezInAmba),ncol(taskGroupfmri)+ncol(networks))
  rownames(resultsTableCorrelation)<-humanEntrezInAmba
  colnames(resultsTableCorrelation)<-c(colnames(taskGroupfmri),colnames(networks))
  
  
  for(k in 1:length(humanEntrezInAmba)){
     resultsTableCorrelation[k,1:ncol(taskGroupfmri)]<-apply(abs(taskGroupfmri),2,function(x){
        cor(x,expressionHuman[,match(humanEntrezInAmba[k],entrezIDsOfRowsHuman)],use ="pairwise.complete.obs",method="spearman")
      })
     
    resultsTableCorrelation[k,(1:ncol(networks))+ncol(taskGroupfmri)]<-apply(abs(networks),2,function(x){
        cor(x,expressionHuman[,match(humanEntrezInAmba[k],entrezIDsOfRowsHuman)],use ="pairwise.complete.obs",method="spearman")
      })
  }
  write.csv2(resultsTableCorrelation,paste0('../storage/',paste0("csvs/",prefix),"gene_task_network_correlation_brainwide.csv"))
 
  head(resultsTableCorrelation)
##          EM-FACES      GA-AVG   GA-REWARD      LA-AVG    MO-RF+LF     MO-RH+LH
## 10043 -0.12225419 -0.13932558 -0.13747772 -0.21584027 -0.09137790 -0.142723895
## 1004  -0.06171345 -0.04170594 -0.03138104 -0.07768244  0.05390198  0.004952909
## 10083 -0.12778216 -0.10381126 -0.10113888 -0.19603992 -0.06891809 -0.115020047
## 10160  0.04232188  0.07811879  0.04979313  0.20333776 -0.07031606 -0.044040218
## 1016   0.19970219  0.18264748  0.18495099  0.26711596  0.16196002  0.221450100
## 10270  0.19452623  0.24256334  0.23319998  0.39914308  0.22099737  0.283286716
##              MO-T       RE-AVG      SO-TOM     WM-REST     WM-FACE
## 10043 -0.09264273 -0.112950052 -0.10688686 -0.16614041 -0.15502639
## 1004   0.01075959 -0.034766388 -0.03058994 -0.02582245 -0.03220408
## 10083 -0.08372677 -0.075708591 -0.15120910 -0.12257321 -0.10260770
## 10160  0.04798651 -0.005703194  0.08324693  0.07407935  0.04643326
## 1016   0.15469942  0.195537701  0.16291591  0.21822855  0.18827356
## 10270  0.20961120  0.188216037  0.16667897  0.29194672  0.26205380
##       cortico-limbic_network mPFC-AMY_network mPFC-Acb_network
## 10043            -0.15934794      -0.06373730     -0.074460375
## 1004              0.11933947      -0.03413417      0.005385308
## 10083            -0.05207671      -0.11742923     -0.093627618
## 10160             0.28324399       0.20443041      0.160412495
## 1016             -0.07207181       0.04023532      0.071913270
## 10270             0.26749063       0.09708940      0.143436205
##       deafult_mode_network salience_network central_executive_network
## 10043          -0.09128415       0.09632715              -0.079311424
## 1004           -0.16226170       0.07869912               0.175555986
## 10083          -0.07432311      -0.05460801              -0.003010283
## 10160           0.29344091       0.12010103               0.010111585
## 1016           -0.01010257      -0.08807999               0.018053726
## 10270           0.16956488      -0.19097613               0.193567471
##       sensorimotor_network visuospatial_network fronto-parietal
## 10043          -0.04617075         -0.025183888     -0.05240495
## 1004            0.09820619          0.006987545     -0.01184755
## 10083           0.00477559         -0.016036088     -0.10640415
## 10160           0.06616453          0.036443746      0.08836849
## 1016            0.10784157          0.115120549      0.11087893
## 10270           0.15316757          0.113855032      0.18731775
##       dorsal_attention ventral_attention
## 10043      -0.05721863      -0.057332473
## 1004       -0.00368351       0.003709607
## 10083      -0.08159040      -0.103653613
## 10160       0.08971130       0.099709028
## 1016        0.10450219       0.103227171
## 10270       0.17103961       0.186488822
  gene_wise_network_correlation<-resultsTableCorrelation
  gene_wise_network_correlationOrg<-gene_wise_network_correlation
  
  gene_wise_network_correlation[,1:11]<-rank(gene_wise_network_correlation[,1:11],ties.method="max",na.last=FALSE)/length(unlist(gene_wise_network_correlation[,1:11]))
  gene_wise_network_correlation[,12:22]<-rank(gene_wise_network_correlation[,12:22],ties.method="max",na.last=FALSE)/length(unlist(gene_wise_network_correlation[,12:22]))
  
  
  gene_wise_network_correlation[,]<-t(apply(gene_wise_network_correlation[,],1,function(x){
   rank(x,ties.method="max",na.last=FALSE)/length(x)
  }))
    
  gene_wise_network_correlation[,]<-t(apply(gene_wise_network_correlation[,],1,function(x){
    x[x==min(x)]<-0
    x
  }))
  
  
  gene_wise_network_correlation[gene_wise_network_correlationOrg<=0.1]<-0
#DNDS to network correlation

  time_clusters_celltype_cor<-cor(dnds_per_gene,(gene_wise_network_correlation),method="spearman",use="pairwise.complete.obs")
    
  rownames(time_clusters_celltype_cor)<-colnames(dnds_per_gene)
  colnames(time_clusters_celltype_cor)<-colnames(gene_wise_network_correlation)
  
  
  time_clusters_celltype_cor<-time_clusters_celltype_cor[,apply(time_clusters_celltype_cor,2,function(x){sum(is.na(x))})<nrow(time_clusters_celltype_cor)]
  time_clusters_celltype_cor[is.na(time_clusters_celltype_cor)]<-0
  
  
  heatmap.2((time_clusters_celltype_cor),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Spearman rank correlation",dendrogram="col",Rowv=NA)

  dev.copy(png,paste0("../storage/biclustering/DNDS_col_Network_col_cor.png"),width=1000,height=500)
## png 
##   3
  dev.off()
## png 
##   2
  write.csv2(time_clusters_celltype_cor,paste0("../storage/biclustering/DNDS_col_Network_col_cor.csv"))
  
  time_clusters_celltype_cor_norm<-t(apply(time_clusters_celltype_cor,1,scale))
  time_clusters_celltype_cor_norm[time_clusters_celltype_cor_norm>3]<-3
  time_clusters_celltype_cor_norm[time_clusters_celltype_cor_norm<(-3)]<-(-3)
  
  colnames(time_clusters_celltype_cor_norm)<-colnames(time_clusters_celltype_cor)
  rownames(time_clusters_celltype_cor_norm)<-rownames(time_clusters_celltype_cor)
  
    heatmap.2((time_clusters_celltype_cor_norm),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Row z-score",dendrogram="col",Rowv=NA)

  dev.copy(png,paste0("../storage/biclustering/DNDS_col_Network_col_cor_row_norm.png"),width=1000,height=500)
## png 
##   3
  dev.off()
## png 
##   2
# Calculate functional maps

############ RUN UNTIL HERE, THEN COMPUTE p-value maps  ############
############ YOU SHOULD BETTER RUN THIS ON A MACHINE    ############
############ WITH >10 CORES SINCE THIS CAN TAKE A WHILE ############

#you don't need to run this, results are included
#source("predict_human_expression_synergy.R")
#Compare functional maps (p-value maps) of branches to taskfMRI

resultsTableCorrelation<-matrix(0,length(transistionsToUseForFunctionalMap),ncol(taskGroupfmri))
rownames(resultsTableCorrelation)<-transistionsToUseForFunctionalMap
colnames(resultsTableCorrelation)<-colnames(taskGroupfmri)

#load pvalue maps an correlate it with teh results table
actRow<-1
for(k in transistionsToUseForFunctionalMap){
  if(file.exists(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds'))){
    pvalsExprHuman<--log10(readRDS(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds')))
    pvalsExprHuman[pvalsExprHuman>3]<-3
    resultsTableCorrelation[actRow,]<-apply(abs(taskGroupfmri),2,function(x){
      cor(x,pvalsExprHuman,use ="pairwise.complete.obs",method="spearman")
    })
    actRow<-actRow+1
  }
}

rownames(resultsTableCorrelation)<-colnames(dnds_per_gene)[transistionsToUseForFunctionalMap]

heatmap.2(t(resultsTableCorrelation),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Spearman rank correlation",dendrogram="row",Colv=NA)

dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/additional_plots/branch_to_task_rank_correlation.png"))
## png 
##   3
dev.off()
## png 
##   2
heatmap.2(apply(resultsTableCorrelation,1,function(x){
    rank(x,ties.method="max",na.last=FALSE)/length(x)
  }),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "white","orange", "red"),key.xlab = "Column rank normalization",dendrogram="row",Colv=NA)

dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/additional_plots/branch_to_task_rank_correlation_column_normalized.png"))
## png 
##   3
dev.off()
## png 
##   2
#Compare branches to networks


resultsTableCorrelation<-matrix(0,length(transistionsToUseForFunctionalMap),ncol(networks))
rownames(resultsTableCorrelation)<-transistionsToUseForFunctionalMap
colnames(resultsTableCorrelation)<-colnames(networks)

actRow<-1
for(k in transistionsToUseForFunctionalMap){
  if(file.exists(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds')) && file.exists(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds'))){
  pvalsExprHuman<--log10(readRDS(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds')))
    pvalsExprHuman[pvalsExprHuman>3]<-3

    resultsTableCorrelation[actRow,]<-apply(abs(networks),2,function(x){
      cor(x,pvalsExprHuman,use ="pairwise.complete.obs",method="spearman")
    })
  }
  actRow<-actRow+1
}
rownames(resultsTableCorrelation)<-colnames(dnds_per_gene)[transistionsToUseForFunctionalMap]


heatmap.2(t(resultsTableCorrelation),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Spearman rank correlation",dendrogram="row",Colv=NA)

dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/additional_plots/branch_to_network_rank_correlation.png"))
## png 
##   3
dev.off()
## png 
##   2
heatmap.2(apply(resultsTableCorrelation,1,function(x){
    rank(x,ties.method="max",na.last=FALSE)/length(x)
  }),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "white","orange", "red"),key.xlab = "Column rank normalization",dendrogram="row",Colv=NA)

dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/additional_plots/branch_to_network_rank_correlation_column_normalized.png"))
## png 
##   3
dev.off()  
## png 
##   2
#Compare branches to networks and tasks combined
resultsTableCorrelation<-matrix(0,length(transistionsToUseForFunctionalMap),ncol(networks)+ncol(taskGroupfmri))
rownames(resultsTableCorrelation)<-transistionsToUseForFunctionalMap
colnames(resultsTableCorrelation)<-c(colnames(networks),colnames(taskGroupfmri))

actRow<-1
for(k in transistionsToUseForFunctionalMap){
  if(file.exists(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds'))){

    pvalsExprHuman<--log10(readRDS(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds')))
    pvalsExprHuman[pvalsExprHuman>3]<-3
    
    resultsTableCorrelation[actRow,]<-c(apply(abs(networks),2,function(x){
      cor(x,pvalsExprHuman,use ="pairwise.complete.obs",method="spearman")
    }),
    apply(abs(taskGroupfmri),2,function(x){
      cor(x,pvalsExprHuman,use ="pairwise.complete.obs",method="spearman")
    })
    )
  }
  actRow<-actRow+1
}

resultsTableCorrelation[,1:ncol(networks)]<-rank(resultsTableCorrelation[,1:ncol(networks)],ties.method="max",na.last=FALSE)/length(unlist(resultsTableCorrelation[,1:ncol(networks)]))
resultsTableCorrelation[,(1:ncol(taskGroupfmri))+ncol(networks)]<-rank(resultsTableCorrelation[,(1:ncol(taskGroupfmri))+ncol(networks)],ties.method="max",na.last=FALSE)/length(unlist(resultsTableCorrelation[,(1:ncol(taskGroupfmri))+ncol(networks)]))
rownames(resultsTableCorrelation)<-colnames(dnds_per_gene)[transistionsToUseForFunctionalMap]

heatmap.2(t(resultsTableCorrelation),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Block-ranks",dendrogram="row",Colv=NA)

dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/branch_to_combined_rank_correlation.png"))
## png 
##   3
dev.off()
## png 
##   2
heatmap.2(t(apply(resultsTableCorrelation,2,function(x){
    rank(x,ties.method="max",na.last=FALSE)/length(x)
  })),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "white","orange", "red"),key.xlab = "Row rank normalization",dendrogram="row",Colv=NA)

dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/additional_plots/branch_to_combined_rank_correlation_row_normalized.png"))  
## png 
##   3
dev.off()
## png 
##   2
heatmap.2(apply(resultsTableCorrelation,1,function(x){
    rank(x,ties.method="max",na.last=FALSE)/length(x)
  }),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "white","orange", "red"),key.xlab = "Column rank normalization",dendrogram="row",Colv=NA)

dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/branch_to_combined_rank_correlation_column_normalized.png"))  
## png 
##   3
dev.off()
## png 
##   2
#Set large brain regions
regionNamesList<-c("myelencephalon","cerebellum","pons","pretectal region","hypothalamus","subthalamus","thalamus","habenular nuclei","paraventricular nuclei of thalamus","pineal gland","basal forebrain","basal ganglia","claustrum","fusiform gyrus","Heschl's gyrus","inferior temporal gyrus","middle temporal gyrus","planum polare","planum temporale","superior temporal gyrus","temporal pole","transverse gyri","inferior parietal lobule","paracentral lobule, posterior part","postcentral gyrus","precuneus","superior parietal lobule","cuneus","inferior occipital gyrus","lingual gyrus","occipital pole","occipito-temporal gyrus","superior occipital gyrus","cingulate gyrus","hippocampal formation","parahippocampal gyrus","piriform cortex","anterior orbital gyrus","frontal operculum","frontal pole","gyrus rectus","inferior frontal gyrus","inferior rostral gyrus","lateral orbital gyrus","medial orbital gyrus","middle frontal gyrus","paracentral lobule, anterior part","parolfactory gyri","posterior orbital gyrus","precentral gyrus","superior frontal gyrus","superior rostral gyrus","insula","amygdalohippocampal transition zone","basolateral nucleus","basomedial nucleus","central nucleus","cortico-medial group","lateral nucleus","central gray substance of midbrain","Edinger-Westphal nucleus","interstitial nucleus of Cajal","midbrain raphe nuclei","midbrain reticular formation","nucleus of Darkschewitsch","oculomotor nuclear complex","red nucleus","substantia nigra","trochlear nucleus","ventral tegmental area","inferior colliculus","superior colliculus")

regionsIDsHumanLarge<-list()
hasBiopsySites<-c()
for(actName in regionNamesList){
  idOfName<-getIDbyName(documentHuman$msg[[1]],actName)$id
  regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-idOfName
  amountOfBiopsySites<-sum(getAtlasRegionsOfIDHuman(atlasRegionsHuman,idOfName)>0)
  if(amountOfBiopsySites==0){
    print(paste(actName,"with id",idOfName,"has no biopsySites!"))
  }
  hasBiopsySites<-c(hasBiopsySites,amountOfBiopsySites>0)
}

print(paste("Length of regionNamesList:",length(regionNamesList)," should be equal to length of regionsIDsHumanLarge:",length(regionsIDsHumanLarge)))
## [1] "Length of regionNamesList: 72  should be equal to length of regionsIDsHumanLarge: 72"
regionNamesList<-regionNamesList[hasBiopsySites]
regionsIDsHumanLarge<-regionsIDsHumanLarge[hasBiopsySites]
#Functions and precomputation for plotting human data

recursiveColorChange<-function(actNode,structureIDs,colorsOfStructures,restWhite=FALSE){
  for(actChild in 1:xmlSize(actNode)){
    attribs<-xmlAttrs(actNode[[actChild]])

    if(!is.na(attribs["style"])){
      isStructure<-attribs["structure_id"]==structureIDs
      if(sum(isStructure,na.rm=TRUE)>0){
       #print(attribs["style"])
        if( !(substr(attribs["style"],nchar(attribs["style"])-6,nchar(attribs["style"]))=="#f2f1f0")){
          if(!(substr(attribs["style"],nchar(attribs["style"])-6,nchar(attribs["style"]))=="#231f20")){
            if(!(substr(attribs["style"],nchar(attribs["style"])-6,nchar(attribs["style"]))=="#6d6e70")){
              if(!substr(attribs["style"],nchar(attribs["style"])-6,nchar(attribs["style"]))=="#872601"){
                if(is.na(colorsOfStructures[isStructure])){
                  if(restWhite) attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
                  else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#cdcdcd")
                }else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),colorsOfStructures[isStructure])
              }else{
                if(restWhite) attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
                else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#231f20")
              }
            }else{
              if(restWhite) attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
              else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#6d6e70")
            }
          }else{
        
            if(restWhite) attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
            else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#231f20")
          }
          
        }else{
          if(restWhite) attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
          else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
        }
        xmlAttrs(actNode[[actChild]])<-attribs
      }else{
        print("BUG")
      }
    }
    if(xmlSize(actNode[[actChild]])>0){
      recursiveColorChange(actNode[[actChild]],structureIDs,colorsOfStructures,restWhite)
    }
    
  }
}

recursiveGetAllStructures<-function(actNode){
  structureIDs<-c()
  for(actChild in 1:xmlSize(actNode)){
    attribs<-xmlAttrs(actNode[[actChild]])
    
    if(!is.na(attribs["style"])){
     structureIDs<-c(attribs["structure_id"],structureIDs)
    }
    if(xmlSize(actNode[[actChild]])>0){
      structureIDs<-unique(c(structureIDs,recursiveGetAllStructures(actNode[[actChild]])))
    }
    
  }
  return(structureIDs)
}

getSmallestParentAtlasRegions <- function(atlasRegionsHuman,ID){
  if(is.null(ID)){
    return(rep(FALSE,length(atlasRegionsHuman)))
  }
  newAtlasRegions<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,ID)
 
  if(sum(newAtlasRegions>0)==0){
    newAtlasRegions<-getSmallestParentAtlasRegions(atlasRegionsHuman,getAcronymByID(documentHuman$msg[[1]],ID)$parent_structure_id)
  }
    
  return(newAtlasRegions)
}

uniqueRegionsHuman<-c()

for (x in 1:104){
  
  fileName <- paste0("../storage/slice_svg/",x,".svg")
  if(!file.exists(fileName)){
    if(!dir.exists("../storage/slice_svg")){
      dir.create("../storage/slice_svg")
    }
    if(!file.exists("../storage/slice_svg/slices.xml")){
     download.file(paste0("http://api.brain-map.org/api/v2/data/query.xml?criteria=model::AtlasImage,rma::criteria,[annotated$eqtrue],alternate_images[image_type$eq%27Atlas+-+Human%27],rma::options[order$eq%27sub_images.section_number%27][num_rows$eqall]"),"../storage/slice_svg/slices.xml")
    }
     
    slicesXML=xmlParse("../storage/slice_svg/slices.xml")
    slicesXMLTop = xmlRoot(slicesXML)
    slicesID<-xmlValue(slicesXMLTop[['atlas-images']][[x]][["id"]][1]$text)
    download.file(paste0("http://api.brain-map.org/api/v2/svg_download/",slicesID,"?groups=265297119,266932194,266932196,266932197"),fileName)
  }
  
  xmlfile=xmlParse(fileName)
  xmltop = xmlRoot(xmlfile)
  
  uniqueRegionsHuman<-unique(c(uniqueRegionsHuman,recursiveGetAllStructures(xmltop[['g']])))
}

uniqueRegionsHumanAtlas<-list()

for(actRegion in 1:length(uniqueRegionsHuman)){
  uniqueRegionsHumanAtlas[[actRegion]]<-getSmallestParentAtlasRegions(atlasRegionsHuman,uniqueRegionsHuman[actRegion])
}

Mode <- function(x) {
  ux <- unique(x[!is.na(x) & x>0])
  if(length(ux)==0){
    return(0)
  }
  ux[which.max(tabulate(match(x[!is.na(x) & x>0], ux)))]
}


plot_mri_style_image_human<-function(path_to_file,plotVals,maxVal,mip,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=FALSE,restWhite=FALSE){

#volume that will be plotted 
m<-array(0, length(uniqueRegionsHuman))

plotVals[plotVals>maxVal]<-maxVal

if(isAlreadyInUniqueRegionsFormat==FALSE){
  for(actRegion in 1:length(uniqueRegionsHuman)){
    if(mip==FALSE)
      m[actRegion]<-Mode(plotVals[uniqueRegionsHumanAtlas[[actRegion]]])
    else m[actRegion]<-mean(plotVals[uniqueRegionsHumanAtlas[[actRegion]]],na.rm=TRUE)
  }
}else{
  m<-plotVals
}

#m[is.nan(m)]<-0
#m[is.na(m)]<-0
if(mip){
  m[m<0]<-0
}else{
  m[m<1]<-NA
}


if(!is.null(path_to_file)){
  png(paste0(path_to_file),width=2000,height=2500)  
}
  par(mfrow=c(4,3),mar = rep(0.1, 4)) 
  
  rb<-c()
  if(discreteColorScale==TRUE){
    rb<-getColorScale(maxVal)
  }else{
    m<-round(m*100)
    rb<-getColorScale(maxVal*100+1)
  }
  
  for (x in 1:11){
    fileName <- paste0("../storage/slice_svg/",x*9,".svg")
    
    xmlfile=xmlParse(fileName)
    xmltop = xmlRoot(xmlfile)
    
    
    if(discreteColorScale==TRUE){
      recursiveColorChange(xmltop[['g']],uniqueRegionsHuman,rb[m],restWhite=restWhite)
    }else{
      recursiveColorChange(xmltop[['g']],uniqueRegionsHuman,rb[m+1],restWhite=restWhite)
    }
    
   
    brainSVG <- rsvg(charToRaw(as(xmlfile, "character")),width=1500)
    plot(c(0,0), type="n", axes=F, xlab="", ylab="")
    rasterImage(brainSVG,1,-1,2,1)
  }
  image.scale <- function(z, zlim, col = heat.colors(12),
                          breaks, horiz=TRUE, ylim=NULL, xlim=NULL, ...){
    if(!missing(breaks)){
      if(length(breaks) != (length(col)+1)){stop("must have one more break than colour")}
    }
    if(missing(breaks) & !missing(zlim)){
      breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1)) 
    }
    if(missing(breaks) & missing(zlim)){
      zlim <- range(z, na.rm=TRUE)
      zlim[2] <- zlim[2]+c(zlim[2]-zlim[1])*(1E-3)#adds a bit to the range in both directions
      zlim[1] <- zlim[1]-c(zlim[2]-zlim[1])*(1E-3)
      breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
    }
    poly <- vector(mode="list", length(col))
    for(i in seq(poly)){
      poly[[i]] <- c(breaks[i], breaks[i+1], breaks[i+1], breaks[i])
    }
    xaxt <- ifelse(horiz, "s", "n")
    yaxt <- ifelse(horiz, "n", "s")
    if(horiz){YLIM<-c(0,1); XLIM<-range(breaks)}
    if(!horiz){YLIM<-range(breaks); XLIM<-c(0,1)}
    if(missing(xlim)) xlim=XLIM
    if(missing(ylim)) ylim=YLIM
    plot(1,1,t="n",ylim=ylim, xlim=xlim, xaxt=xaxt, yaxt=yaxt, xaxs="i", yaxs="i",cex.axis=6, ...)  
    for(i in seq(poly)){
      if(horiz){
        polygon(poly[[i]], c(0,0,1,1), col=col[i], border=NA)
      }
      if(!horiz){
        polygon(c(0,0,1,1), poly[[i]], col=col[i], border=NA)
      }
    }

  }
  
  par(mgp=c(5,6,0))
  if(!is.null(path_to_file)){
    par(mar=c(15,5,15,5))
  }else{
    par(mar=c(3,1,3,1))
  }

   if(mip){
     if(discreteColorScale==TRUE){
       image.scale(c(0,((maxVal))),col=getColorScale((maxVal)))
     }else{
       image.scale(c(0,((maxVal))),col=getColorScale((maxVal*100)))
     }
      
    }else{
     if(discreteColorScale==TRUE){
       image.scale(c(1,((maxVal))),col=getColorScale((maxVal)))
     }else{
       image.scale(c(1,((maxVal))),col=getColorScale((maxVal*100)))
     }
    }
  
  
  if(!is.null(path_to_file)){
    dev.off()
  }

}

plot_mri_style_image_human_all_slices<-function(path_to_file,plotVals,maxVal,mip,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=FALSE,restWhite=FALSE){

#volume that will be plotted 
m<-array(0, length(uniqueRegionsHuman))

plotVals[plotVals>maxVal]<-maxVal

if(isAlreadyInUniqueRegionsFormat==FALSE){
  for(actRegion in 1:length(uniqueRegionsHuman)){
    if(mip==FALSE)
      m[actRegion]<-Mode(plotVals[uniqueRegionsHumanAtlas[[actRegion]]])
    else m[actRegion]<-mean(plotVals[uniqueRegionsHumanAtlas[[actRegion]]],na.rm=TRUE)
  }
}else{
  m<-plotVals
}

#m[is.nan(m)]<-0
#m[is.na(m)]<-0
if(mip){
  m[m<0]<-0
}else{
  m[m<1]<-NA
}


  png(paste0(path_to_file),width=10000,height=4000, type="cairo")  
  par(mfrow=c(6,17),mar = rep(0.1, 4)) 
  
  rb<-c()
  if(discreteColorScale==TRUE){
    rb<-getColorScale(maxVal)
  }else{
    m<-round(m*100)
    rb<-getColorScale(maxVal*100+1)
  }
  
   for (x in 1:101){
      fileName <- paste0("../storage/slice_svg/",x,".svg")
    
      
      xmlfile=xmlParse(fileName)
      xmltop = xmlRoot(xmlfile)
      
      
      if(discreteColorScale==TRUE){
        recursiveColorChange(xmltop[['g']],uniqueRegionsHuman,rb[m],restWhite=restWhite)
      }else{
        recursiveColorChange(xmltop[['g']],uniqueRegionsHuman,rb[m+1],restWhite=restWhite)
      }
      
      brainSVG <- rsvg(charToRaw(as(xmlfile, "character")),width=1000)
      plot(c(0,0), type="n", axes=F, xlab="", ylab="")
      rasterImage(brainSVG,1,-1,2,1)
    }
  image.scale <- function(z, zlim, col = heat.colors(12),
                          breaks, horiz=TRUE, ylim=NULL, xlim=NULL, ...){
    if(!missing(breaks)){
      if(length(breaks) != (length(col)+1)){stop("must have one more break than colour")}
    }
    if(missing(breaks) & !missing(zlim)){
      breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1)) 
    }
    if(missing(breaks) & missing(zlim)){
      zlim <- range(z, na.rm=TRUE)
      zlim[2] <- zlim[2]+c(zlim[2]-zlim[1])*(1E-3)#adds a bit to the range in both directions
      zlim[1] <- zlim[1]-c(zlim[2]-zlim[1])*(1E-3)
      breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
    }
    poly <- vector(mode="list", length(col))
    for(i in seq(poly)){
      poly[[i]] <- c(breaks[i], breaks[i+1], breaks[i+1], breaks[i])
    }
    xaxt <- ifelse(horiz, "s", "n")
    yaxt <- ifelse(horiz, "n", "s")
    if(horiz){YLIM<-c(0,1); XLIM<-range(breaks)}
    if(!horiz){YLIM<-range(breaks); XLIM<-c(0,1)}
    if(missing(xlim)) xlim=XLIM
    if(missing(ylim)) ylim=YLIM
    plot(1,1,t="n",ylim=ylim, xlim=xlim, xaxt=xaxt, yaxt=yaxt, xaxs="i", yaxs="i",cex.axis=6, ...)  
    for(i in seq(poly)){
      if(horiz){
        polygon(poly[[i]], c(0,0,1,1), col=col[i], border=NA)
      }
      if(!horiz){
        polygon(c(0,0,1,1), poly[[i]], col=col[i], border=NA)
      }
    }

  }
  
   par(mgp=c(3,6,0))
   par(mar=c(15,5,15,5))

   if(mip){
     if(discreteColorScale==TRUE){
       image.scale(c(0,((maxVal))),col=getColorScale((maxVal)))
     }else{
       image.scale(c(0,((maxVal))),col=getColorScale((maxVal*100)))
     }
      
    }else{
     if(discreteColorScale==TRUE){
       image.scale(c(1,((maxVal))),col=getColorScale((maxVal)))
     }else{
       image.scale(c(1,((maxVal))),col=getColorScale((maxVal*100)))
     }
    }
  
   dev.off()

}

plot_mri_style_image_human_all_slices_and_significance<-function(path_to_file,plotPvals,fdr005,fdr01,justDrawOriginalColors=FALSE){
    getColorScale<-function(amountOfColors,threshold){
        greyColor<-colorpanel(amountOfColors,"#f9f9f9","#a8a8a8","#7a7a7a")
        heatColor<-colorpanel(amountOfColors,"#ffffb2","#feb24c","#bd0026")
        
        if(is.na(threshold))
          return(greyColor)
        if(threshold<=1 || threshold==Inf)
          return(greyColor)
        return(c(greyColor[1:floor(threshold)],heatColor[(floor(threshold)+1):amountOfColors]))
      }
  
    #volume that will be plotted 
    m<-array(0, length(uniqueRegionsHuman))
    plotPvals[is.na(plotPvals)]<-1
    
    
    for(actRegion in 1:length(uniqueRegionsHuman)){
      m[actRegion]<-quantile(-log10(plotPvals[uniqueRegionsHumanAtlas[[actRegion]]]),0.9,na.rm=TRUE)
    }
    
    m[is.nan(m)]<-0
    m[is.na(m)]<-0
    m[m<0]<-0
    
    minimumPval<--log10(10^(-8))
    
    m[m>minimumPval]<-minimumPval
    
    png(paste0(path_to_file),width=10000,height=4000, type="cairo")  
    par(mfrow=c(6,17),mar = rep(0.1, 4)) 
    
    
    m<-round(m*1000)
    rb<-getColorScale(minimumPval*1000+1,-log10(fdr01)*1000+1)
    
    for (x in 1:101){
      fileName <- paste0("../storage/slice_svg/",x,".svg")
    
      
      xmlfile=xmlParse(fileName)
      xmltop = xmlRoot(xmlfile)
      if(justDrawOriginalColors==FALSE){
        recursiveColorChange(xmltop[['g']],uniqueRegionsHuman,rb[m+1])
      }
      
      brainSVG <- rsvg(charToRaw(as(xmlfile, "character")),width=1000)
      plot(c(0,0), type="n", axes=F, xlab="", ylab="")
      rasterImage(brainSVG,1,-1,2,1)
    }
    image.scale <- function(z, zlim, col = heat.colors(12),
                            breaks, horiz=TRUE, ylim=NULL, xlim=NULL, ...){
      if(!missing(breaks)){
        if(length(breaks) != (length(col)+1)){stop("must have one more break than colour")}
      }
      if(missing(breaks) & !missing(zlim)){
        breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1)) 
      }
      if(missing(breaks) & missing(zlim)){
        zlim <- range(z, na.rm=TRUE)
        zlim[2] <- zlim[2]+c(zlim[2]-zlim[1])*(1E-3)#adds a bit to the range in both directions
        zlim[1] <- zlim[1]-c(zlim[2]-zlim[1])*(1E-3)
        breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
      }
      poly <- vector(mode="list", length(col))
      for(i in seq(poly)){
        poly[[i]] <- c(breaks[i], breaks[i+1], breaks[i+1], breaks[i])
      }
      xaxt <- ifelse(horiz, "s", "n")
      yaxt <- ifelse(horiz, "n", "s")
      if(horiz){YLIM<-c(0,1); XLIM<-range(breaks)}
      if(!horiz){YLIM<-range(breaks); XLIM<-c(0,1)}
      if(missing(xlim)) xlim=XLIM
      if(missing(ylim)) ylim=YLIM
      plot(1,1,t="n",ylim=ylim, xlim=xlim, xaxt=xaxt, yaxt=yaxt, xaxs="i", yaxs="i",cex.axis=8, ...)  
      for(i in seq(poly)){
        if(horiz){
          polygon(poly[[i]], c(0,0,1,1), col=col[i], border=NA)
        }
        if(!horiz){
          polygon(c(0,0,1,1), poly[[i]], col=col[i], border=NA)
        }
      }
      if(!(fdr005==0)){
        lines(c(-log10(fdr005),-log10(fdr005)),c(0,1),lty=3,lwd=5)
      }
      if(!(fdr01==0)){
        lines(c(-log10(fdr01),-log10(fdr01)),c(0,1),lwd=5)
      }
    }
    
    par(mgp=c(3,6,0))
    par(mar=c(15,5,15,5))
    
    image.scale(c(0,minimumPval),col=rb)
    dev.off()
    
  }
#Calcualte evolutionary age of every biopsy site in human brain

pvalueMapsRessource<-matrix(0, length(transistionsToUseForFunctionalMap),sum(atlasRegionsHuman>0))
humanFunctionalMapsLogPvalueTop<-matrix(0, length(transistionsToUseForFunctionalMap),sum(atlasRegionsHuman>0))
humanFunctionalMapsUniqueMapIsSignificant<-matrix(FALSE, length(transistionsToUseForFunctionalMap),length(uniqueRegionsHumanAtlas))

 plot_mri_style_image_human_all_slices_and_significance(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_region_colors.png"),rep(0,ncol(humanFunctionalMapsLogPvalueTop)),0,0,justDrawOriginalColors=TRUE)
## png 
##   2
branch<-1
for(k in transistionsToUseForFunctionalMap){
  print(paste0("branch: ",branch))

  humanFunctionalMapsLogPvalueTop[branch,]<-readRDS(paste0(functionalMapsHuman,'/','pvalues/',"dndsweighted_trans",k,'_ttest_geneExpr.rds'))
  
  sigThresh<-(-log10(max(humanFunctionalMapsLogPvalueTop[branch,][p.adjust(humanFunctionalMapsLogPvalueTop[branch,],method="BH")<=0.1])))
  for(actRegion in 1:length(uniqueRegionsHumanAtlas)){
    if(!is.na(sigThresh)){
      if(quantile(-log10(humanFunctionalMapsLogPvalueTop[branch,][uniqueRegionsHumanAtlas[[actRegion]]]),0.9,na.rm=TRUE)>=sigThresh){
        humanFunctionalMapsUniqueMapIsSignificant[branch,actRegion]<-TRUE
      }
    }
  }
  
  fdr005<-max(humanFunctionalMapsLogPvalueTop[branch,][p.adjust(humanFunctionalMapsLogPvalueTop[branch,],method="BH")<=0.05])
  fdr01<-max(humanFunctionalMapsLogPvalueTop[branch,][p.adjust(humanFunctionalMapsLogPvalueTop[branch,],method="BH")<=0.1])
  plot_mri_style_image_human_all_slices_and_significance(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_",colsToUseNames[k],"_90percentile_significant_sites",transistionsToUseForFunctionalMapPostfix,".png"),humanFunctionalMapsLogPvalueTop[branch,],fdr005,fdr01)
  
  pvalueMapsRessource[branch,]<-humanFunctionalMapsLogPvalueTop[branch,]
  humanFunctionalMapsLogPvalueTop[branch,]<-(-log10(humanFunctionalMapsLogPvalueTop[branch,]))
  branch<-branch+1
}
## [1] "branch: 1"
## [1] "branch: 2"
## [1] "branch: 3"
## [1] "branch: 4"
## [1] "branch: 5"
## [1] "branch: 6"
## [1] "branch: 7"
## [1] "branch: 8"
evolutionaryAgeUniqueRegionsHumanAtlasWiseTop<-unlist(sapply(uniqueRegionsHumanAtlas,function(x){
  if(sum(x)>1){
    which.max(apply(humanFunctionalMapsLogPvalueTop[,x],1,function(y){
      mean(y,na.rm=TRUE)
    }))
  }else{
    which.max(humanFunctionalMapsLogPvalueTop[,x])
  }
}))


humanFunctionalMapsLogPvalueTop[humanFunctionalMapsLogPvalueTop>15]<-15


pvalueMapsRessource<-t(pvalueMapsRessource)
colnames(pvalueMapsRessource)<-colsToUseNames[transistionsToUseForFunctionalMap]

appendSheet<-FALSE
options(scipen = 999)
options(warn=-1)
if(file.exists( paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/pvalue_ressource",transistionsToUseForFunctionalMapPostfix,".xlsx"))){
  file.remove( paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/pvalue_ressource",transistionsToUseForFunctionalMapPostfix,".xlsx"))
}
## [1] TRUE
for(i in 1:ncol(pvalueMapsRessource)){
  
  regionPvals<-cbind(regionInfo,pvalueMapsRessource[,i])
  colnames(regionPvals)[ncol(regionPvals)]<-"pvalue"
  write.xlsx(regionPvals, file = paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/pvalue_ressource",transistionsToUseForFunctionalMapPostfix,".xlsx"),
      sheetName = colnames(pvalueMapsRessource)[i], append = appendSheet,row.names = FALSE)
  
    appendSheet<-TRUE
}
options(warn=0)



evolutionaryAgeBiopsyWiseTop<-unlist(apply((humanFunctionalMapsLogPvalueTop),2,function(x){return(which.max(x))}))
evolutionaryAgeBiopsyWiseTop<-sapply(uniqueRegionsHumanAtlas,function(x){Mode(evolutionaryAgeBiopsyWiseTop[x])})
evolutionaryAgeBiopsyWiseTopSign<-sapply(1:length(evolutionaryAgeBiopsyWiseTop),function(x){
  if(is.na(evolutionaryAgeBiopsyWiseTop[x])){
    return(NA)
  }
  if(humanFunctionalMapsUniqueMapIsSignificant[evolutionaryAgeBiopsyWiseTop[x],x]==FALSE){
    return(NA)
  }else{
    return(evolutionaryAgeBiopsyWiseTop[x])
  }
})


evolutionaryAgeBiopsyWiseTopHom<-unlist(apply((humanFunctionalMapsLogPvalueTop[(length(transistionsToUseForFunctionalMap)-2):length(transistionsToUseForFunctionalMap),]),2,function(x){return(which.max(x)+(length(transistionsToUseForFunctionalMap)-3))}))
evolutionaryAgeBiopsyWiseTopHom<-sapply(uniqueRegionsHumanAtlas,function(x){Mode(evolutionaryAgeBiopsyWiseTopHom[x])})
evolutionaryAgeBiopsyWiseTopHomSign<-sapply(1:length(evolutionaryAgeBiopsyWiseTopHom),function(x){
  if(is.na(evolutionaryAgeBiopsyWiseTopHom[x])){
    return(NA)
  }
  if(humanFunctionalMapsUniqueMapIsSignificant[evolutionaryAgeBiopsyWiseTopHom[x],x]==FALSE){
    return(NA)
  }else{
    return(evolutionaryAgeBiopsyWiseTopHom[x])
  }
})


evolutionaryAgeBiopsyWiseWATop<-apply((humanFunctionalMapsLogPvalueTop),2,function(x){
  trans<-1:length(x)
  if(sum(is.na(x))==length(x)){
    return(NA)
  }
  y<-x
  y[y<0]<-0
  
  return(sum((trans[!is.na(y)])*y[!is.na(y)]/(sum(y,na.rm=TRUE))))
})
evolutionaryAgeBiopsyWiseWATop<-sapply(uniqueRegionsHumanAtlas,function(x){mean(evolutionaryAgeBiopsyWiseWATop[x],na.rm=TRUE)})
#Cluster Evo History By Gene Expression

dir.create(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted"))
## Warning in dir.create(paste0("../storage/",
## prefix, "evoHistoryHuman_dndsweighted")): '..
## \storage\NGT_SPLIT_TREE_9000_NAto0_evoHistoryHuman_dndsweighted' already exists
expressionHumanTopLargeRegions<-sapply(1:length(regionsIDsHumanLarge),function(actRegion){
  actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[[actRegion]][1])>0
  if(sum(actRegionAtlas)==1){
    return(expressionHuman[actRegionAtlas,is.element(entrezIDsOfRowsHuman,humanEntrezInAmba)])
  }else{
    return(apply(expressionHuman[actRegionAtlas,is.element(entrezIDsOfRowsHuman,humanEntrezInAmba)],2,function(x){
      mean(x,na.rm=TRUE)
    }))
  }
})


colnames(expressionHumanTopLargeRegions)<-1:ncol(expressionHumanTopLargeRegions)
for(actRegion in 1:length(regionsIDsHumanLarge)){
  colnames(expressionHumanTopLargeRegions)[actRegion]<-getAcronymByID(documentHuman$msg[[1]],regionsIDsHumanLarge[[actRegion]])$acronym
}

 cordist<-as.dist((1-cor(expressionHumanTopLargeRegions))/2)
 cordist[is.na(cordist)]<-0
 hc <- hclust(cordist,method="ward.D2")

#hc <- hclust(dist(t(expressionHumanTopLargeRegions)),method="ward.D2")
clustersExpression<-cutree(hc, k = 8)
hist(clustersExpression)

plot(hc)

dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/gene_expression_clustering_dendrogram",transistionsToUseForFunctionalMapPostfix,".png"),width=900,height=400)
## png 
##   3
dev.off()
## png 
##   2
clustersExpressionCSV<-paste("Cluster",clustersExpression)
names(clustersExpressionCSV)<-colnames(expressionHumanTopLargeRegions)
sortedclustersExpressionCSV<-sort(clustersExpressionCSV)
write.csv2(sortedclustersExpressionCSV,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/expression_regions_cluster",transistionsToUseForFunctionalMapPostfix,".csv"))

clustersExpressionBiopsyLevel<-rep(NA,length(atlasRegionsHuman))
for(actRegion in 1:length(regionsIDsHumanLarge)){
  actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[[actRegion]][1])
  clustersExpressionBiopsyLevel[actRegionAtlas]<-clustersExpression[actRegion]
}

clusterProfiles<-matrix(0,dim(expressionHumanTopLargeRegions)[1],max(clustersExpression,na.rm=TRUE))
for(timepoint in 1:dim(expressionHumanTopLargeRegions)[1]){
  for(cluster in unique(clustersExpression)){
    clusterProfiles[timepoint,cluster]<-mean(expressionHumanTopLargeRegions[timepoint,clustersExpression==cluster],rm.na=TRUE)
  }
}
library(RColorBrewer)
getColorScale<-function(amountOfColors){
    #return(gsub('.{2}$', '',rainbow(amountOfColors)))
  return(brewer.pal(n = 8, name = "Set2"))
}
rb<-getColorScale(max(clustersExpression,na.rm=TRUE))

plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/expression_clustering",transistionsToUseForFunctionalMapPostfix,".png"),clustersExpressionBiopsyLevel,max(clustersExpressionBiopsyLevel,na.rm=TRUE),mip=FALSE,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=TRUE,restWhite=FALSE)
## png 
##   2
plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/expression_clustering",transistionsToUseForFunctionalMapPostfix,".png"),clustersExpressionBiopsyLevel,max(clustersExpressionBiopsyLevel,na.rm=TRUE),mip=FALSE,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=TRUE,restWhite=FALSE)
## png 
##   2
#Cluster Evo History By Functional Maps

humanFunctionalMapsLogPvalueTopLargeRegions<-sapply(1:length(regionsIDsHumanLarge),function(actRegion){
  actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[[actRegion]][1])>0
  if(sum(actRegionAtlas)==1){
    return(humanFunctionalMapsLogPvalueTop[,actRegionAtlas])
  }else{
    return(apply(humanFunctionalMapsLogPvalueTop[,actRegionAtlas],1,function(x){
      mean(x,na.rm=TRUE)
    }))
  }
})


humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA<-humanFunctionalMapsLogPvalueTopLargeRegions
humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA[is.na(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)]<-0
humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA<-t(apply(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA,1,rank)-1)
humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA<-humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA/apply(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA,1,max)
humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA[humanFunctionalMapsLogPvalueTopLargeRegions==0 | is.na(humanFunctionalMapsLogPvalueTopLargeRegions)]<-0

colnames(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)<-1:ncol(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)
for(actRegion in 1:length(regionsIDsHumanLarge)){
  colnames(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)[actRegion]<-getAcronymByID(documentHuman$msg[[1]],regionsIDsHumanLarge[[actRegion]])$acronym
}

write.csv2(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_region_pvals",transistionsToUseForFunctionalMapPostfix,".csv"))

cordist<-as.dist((1-cor(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA))/2)
cordist[is.na(cordist)]<-0
hc <- hclust(dist(t(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)),method="ward.D2")
clustersEvoHist<-cutree(hc, k = 8)
hist(clustersEvoHist)

plot(hc)

dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_dendrogram",transistionsToUseForFunctionalMapPostfix,".png"),width=900,height=400)
## png 
##   3
dev.off()
## png 
##   2
clustersEvoHistCSV<-paste("Cluster",clustersEvoHist)
names(clustersEvoHistCSV)<-colnames(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)
sortedclustersExpressionCSV<-sort(clustersEvoHistCSV)
write.csv2(sortedclustersExpressionCSV,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_cluster",transistionsToUseForFunctionalMapPostfix,".csv"))

clustersEvoHistBipsyLevel<-rep(NA,length(atlasRegionsHuman))
for(actRegion in 1:length(regionsIDsHumanLarge)){
  actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[[actRegion]][1])
  clustersEvoHistBipsyLevel[actRegionAtlas]<-clustersEvoHist[actRegion]
}

clusterProfiles<-matrix(0,dim(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)[1],max(clustersEvoHist,na.rm=TRUE))
for(timepoint in 1:dim(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)[1]){
  for(cluster in unique(clustersEvoHist)){
    clusterProfiles[timepoint,cluster]<-mean(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA[timepoint,clustersEvoHist==cluster],rm.na=TRUE)
  }
}
library(RColorBrewer)
getColorScale<-function(amountOfColors){
    #return(gsub('.{2}$', '',rainbow(amountOfColors)))
  return(brewer.pal(n = 8, name = "Set2"))
}
rb<-getColorScale(max(clustersEvoHist,na.rm=TRUE))

matplot(1:nrow(clusterProfiles), clusterProfiles, lty=1,type='l', xlab='DNDS column', ylab='rank normalized mean p-value (selection pressure)', col=rb[1:max(clustersEvoHist,na.rm=TRUE)],ylim=c(0,1.1),lwd=2,xaxt="n")
axis(side=1,at=1:nrow(clusterProfiles),labels=transistionsToUseForFunctionalMap)

rownames(clusterProfiles)<-transistionsToUseForFunctionalMap
colnames(clusterProfiles)<-paste("Cl",1:ncol(clusterProfiles))
write.csv2(clusterProfiles,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_timeline",transistionsToUseForFunctionalMapPostfix,".csv"))

legend('topright', inset=.05, legend=paste("Cl",1:ncol(clusterProfiles)), 
                             lty=1, horiz=TRUE, col=rb[1:max(clustersEvoHist,na.rm=TRUE)])

dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_timeline",transistionsToUseForFunctionalMapPostfix,".png"),width=900,height=400)
## png 
##   3
dev.off()
## png 
##   2
plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering",transistionsToUseForFunctionalMapPostfix,".png"),clustersEvoHistBipsyLevel,max(clustersEvoHistBipsyLevel,na.rm=TRUE),mip=FALSE,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=TRUE,restWhite=FALSE)
## png 
##   2
plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_clustering",transistionsToUseForFunctionalMapPostfix,".png"),clustersEvoHistBipsyLevel,max(clustersEvoHistBipsyLevel,na.rm=TRUE),mip=FALSE,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=TRUE,restWhite=FALSE)
## png 
##   2
regionCorrelation<-cor(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)

heatmap.2(regionCorrelation,scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Correlation",Rowv=order.dendrogram(as.dendrogram(hc)),Colv=order.dendrogram(as.dendrogram(hc)),hclustfun=function(x){hc})

dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_correlation",transistionsToUseForFunctionalMapPostfix,".png"))  
## png 
##   3
dev.off()
## png 
##   2
heatmap.2(cor(expressionHumanTopLargeRegions),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Correlation",Rowv=order.dendrogram(as.dendrogram(hc)),Colv=order.dendrogram(as.dendrogram(hc)),hclustfun=function(x){hc})

dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_expression_correlation",transistionsToUseForFunctionalMapPostfix,".png"))  
## png 
##   3
dev.off()
## png 
##   2
heatmap.2(regionCorrelation-cor(expressionHumanTopLargeRegions),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Correlation Distance",Rowv=order.dendrogram(as.dendrogram(hc)),Colv=order.dendrogram(as.dendrogram(hc)),hclustfun=function(x){hc})

dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_expression_correlation_dist",transistionsToUseForFunctionalMapPostfix,".png"))  
## png 
##   3
dev.off()
## png 
##   2
#Network Figure

colorsOfRegions<-sapply(1:length(regionsIDsHumanLarge),function(actRegion){
  paste0("#",getAcronymByID(documentHuman$msg[[1]],regionsIDsHumanLarge[[actRegion]][1])$color_hex_triplet)
})


regionCorrelation<-cor(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)
write.csv2(regionCorrelation,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_correlation.csv"))
links <- matrix(0,nrow=0,ncol=3)


for(i in 1:nrow(regionCorrelation)){
  for(j in i:ncol(regionCorrelation)){
    if(regionCorrelation[i,j]>0 && i!=j){
      links <- rbind(links,c(rownames(regionCorrelation)[i],rownames(regionCorrelation)[j],regionCorrelation[i,j]))
    }
  }
}

links<-data.frame(
    source=links[,1],
    target=links[,2],
    importance=as.numeric(links[,3])*100,
    stringsAsFactors=FALSE
    )

nodes <- data.frame(
    name=rownames(regionCorrelation)
    )
 
# Turn it into igraph object
network <- graph_from_data_frame(d=links, vertices=nodes, directed=F) 
E(network)$color<-gray.colors(100,start = 1, end = 0.1)[round(links$importance)+1]
E(network)$weights<-links$importance
V(network)$frame.color <- "white"
V(network)$label.color <- "black"
V(network)$label.font=2
V(network)$label.family="sans"

layoutingAlg<-layout_with_kk(network)

# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_fr(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_gem(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_graphopt(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_kk(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_lgl(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_mds(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_dh(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_sugiyama(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_nicely(network))


plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layoutingAlg)

plot(network, vertex.color=rb[clustersEvoHist],vertex.size=8, layout=layoutingAlg)

nodecolors<-cbind(rownames(regionCorrelation),colorsOfRegions,rb[clustersEvoHist])
colnames(nodecolors)<-c("regionname","AHBA color","clustercolor")
write.csv2(nodecolors,"C:\\Users\\ganglberger\\Desktop\\nodecolors.csv")

# Make the plot

svg(filename=paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_correlation_graph",transistionsToUseForFunctionalMapPostfix,".svg"),width=20, height=20)
plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layoutingAlg)
dev.off()
## png 
##   2
svg(filename=paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_correlation_graph_cluster_cols",transistionsToUseForFunctionalMapPostfix,".svg"),width=20, height=20)
plot(network, vertex.color=rb[clustersEvoHist],vertex.size=8, layout=layoutingAlg)
dev.off()
## png 
##   2
networksCombined<-cbind(taskGroupfmri,networks)
colnames(networksCombined)<-c(colnames(taskGroupfmri),colnames(networks))

networksCombined_on_network_res<-sapply(1:length(regionsIDsHumanLarge),function(actRegion){
  actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[[actRegion]][1])>0
  if(sum(actRegionAtlas)==1){
    return(networksCombined[actRegionAtlas,])
  }else{
    return(apply(networksCombined[actRegionAtlas,],2,function(x){
      max(x,na.rm=TRUE)
    }))
  }
})

networksCombined_on_network_res[networksCombined_on_network_res<0]<-0
networkscombined_on_network_res_rank<-t(apply(networksCombined_on_network_res,1,function(x){rank(x,ties.method="max",na.last="keep")/length(x)}))
networkscombined_on_network_res_rank[networksCombined_on_network_res==0]<-0

for(act_network_index in 1:ncol(networksCombined)){
  svg(filename=paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_correlation_graph_",colnames(networksCombined)[act_network_index],"",transistionsToUseForFunctionalMapPostfix,".svg"),width=20, height=20)
  
  newColors<-rb[clustersEvoHist]
  newColors[networkscombined_on_network_res_rank[act_network_index,]==0]<-"#FFFFFF"
  
plot(network, vertex.color=newColors,vertex.size=8, layout=layoutingAlg)
dev.off()
  
}
#Plot evolutionary age of every biopsy site in human brain

number_colnames<-1:length(transistionsToUseForFunctionalMap)
names(number_colnames)<-colnames(dnds_per_gene)[transistionsToUseForFunctionalMap]
write.csv2(number_colnames,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_numbers_to_columnnames",transistionsToUseForFunctionalMapPostfix,".csv"))


  getColorScale<-function(amountOfColors){
    return(jet.colors(amountOfColors))
  }
  
    plot_mri_style_image_human(NULL,evolutionaryAgeUniqueRegionsHumanAtlasWiseTop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)

    plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_regionwise_peak",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeUniqueRegionsHumanAtlasWiseTop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png 
##   2
  plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_peak_modeRegion",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseTop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png 
##   2
   plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_peak_modeRegion",transistionsToUseForFunctionalMapPostfix,"_sign.png"),evolutionaryAgeBiopsyWiseTopSign,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png 
##   2
  plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistoryWA_meanRegion",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseWATop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=FALSE,isAlreadyInUniqueRegionsFormat=TRUE)
## png 
##   2
      plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_regionwise_peak",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeUniqueRegionsHumanAtlasWiseTop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png 
##   2
    plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_peak_modeRegion",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseTop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png 
##   2
   plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_peak_modeRegion",transistionsToUseForFunctionalMapPostfix,"_sign.png"),evolutionaryAgeBiopsyWiseTopSign,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png 
##   2
  plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistoryWA_meanRegion",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseWATop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=FALSE,isAlreadyInUniqueRegionsFormat=TRUE)
## png 
##   2
  getColorScale<-function(amountOfColors){
    return(colorpanel(amountOfColors,"#00007F","#F4F4F4","#7F0000"))
  }

  plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistoryWA_meanRegion_br",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseWATop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=FALSE,isAlreadyInUniqueRegionsFormat=TRUE)
## png 
##   2
    plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistoryWA_meanRegion_br",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseWATop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=FALSE,isAlreadyInUniqueRegionsFormat=TRUE)
## png 
##   2
  getColorScale<-function(amountOfColors){
    return(colorpanel(amountOfColors,"#ffffb2","#feb24c","#bd0026"))
  }

  # maxVal<-quantile(humanFunctionalMapsLogPvalueTop[humanFunctionalMapsLogPvalueTop<Inf],0.95,na.rm=TRUE)
  # for(actTrans in 1:length(transistionsToUseForFunctionalMap)){
  #    plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_t",transistionsToUseForFunctionalMap[actTrans],"_mean_of_significant_sites.png"),humanFunctionalMapsLogPvalueTop[actTrans,],maxVal,mip=TRUE)
  # }
#Calculate biclustering

############ RUN UNTIL HERE, THEN COMPUTE BICLUSTERING  ############
############ YOU SHOULD BETTER RUN THIS ON A MACHINE    ############
############ WITH >10 CORES SINCE THIS CAN TAKE A WHILE ############

#you don't need to run this, results are included
#source("GABI_BICLUSTERING.R")
for(actBiclustering in c("ALL_DNDS_","C_mouse_to_human_","E2_hominids_speciation_with_ancestorsNOG_","C3_mouse_to_human_late_no_gorgo_")){
  for(lowOrHigh in c("high","low")){
  
    if(lowOrHigh=="high"){
      biclustering_results<-readRDS(paste0("../storage/biclustering_results/",actBiclustering,"ALL_NETWORKS_0_9_0_75.rds"))
    }else{
      biclustering_results<-readRDS(paste0("../storage/biclustering_results/",actBiclustering,"ALL_NETWORKS_0_1_0_75.rds"))
    }
    if(actBiclustering=="ALL_DNDS_"){
      time_to_use <- c(1:10,11, 12, 13, 14, 15, 18, 19, 16, 17, 20, 21, 10)
    
    #remove unstable clusters
    if(lowOrHigh=="high"){
    }else{
       biclustering_results<-biclustering_results[-c(21)]
    }
    }
    
    if(actBiclustering=="C_mouse_to_human_"){
      time_to_use <- c(11, 12, 13, 14, 15, 18, 19, 16, 17, 20, 21, 10)
    
      #remove unstable clusters
      if(lowOrHigh=="high"){
      }else{
        biclustering_results<-biclustering_results[-c(21)]
      }
    }
    
    if(actBiclustering=="E2_hominids_speciation_with_ancestorsNOG_"){
      time_to_use <- c(18, 15, 7, 19, 16, 17, 20, 21, 8, 9, 10)
      
      #remove unstable clusters
      if(lowOrHigh=="high"){
         biclustering_results<-biclustering_results[-c(6,10)]
      }else{
        biclustering_results<-biclustering_results[-c(21)]
      }
    
    }
    
    if(actBiclustering=="C3_mouse_to_human_late_no_gorgo_"){
    time_to_use <- c(10, 16, 17, 20, 21)
    
    #remove unstable clusters
    if(lowOrHigh=="high"){
      }else{
        biclustering_results<-biclustering_results[-c(15)]
      }
    }  
    
    biclusteringSummary <- c()
    clusterFeatures <- lapply(biclustering_results, function(x) {
    dndsCols <- x$features[x$features <= ncol(normalized_dnds_per_gene)]
    networkCols <- x$features[x$features > ncol(normalized_dnds_per_gene)]
    c(
      paste(colnames(normalized_dnds_per_gene)[dndsCols], collapse = " "),
      paste(colnames(gene_wise_network_correlation)[networkCols - ncol(normalized_dnds_per_gene)], collapse =
              " "),
      length(x$samples),
      paste(apply(as.matrix(normalized_dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
        round(mean(x, na.rm = TRUE), digits = 3)
      }), collapse = " "),
      paste(apply(as.matrix(normalized_dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
        round(sd(x, na.rm = TRUE), digits = 3)
      }), collapse = " "),
      paste(apply(as.matrix(normalized_dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
        round(quantile(x, na.rm = TRUE, probs = 0.1), digits = 3)
      }), collapse = " "),
      paste(apply(as.matrix(normalized_dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
        round(quantile(x, na.rm = TRUE, probs = 0.9), digits = 3)
      }), collapse = " "),
      paste(apply(as.matrix(dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
        round(mean(x, na.rm = TRUE), digits = 3)
      }), collapse = " "),
      paste(apply(as.matrix(dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
        round(sd(x, na.rm = TRUE), digits = 3)
      }), collapse = " "),
      paste(apply(as.matrix(dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
        round(quantile(x, na.rm = TRUE, probs = 0.1), digits = 3)
      }), collapse = " "),
      paste(apply(as.matrix(dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
        round(quantile(x, na.rm = TRUE, probs = 0.9), digits = 3)
      }), collapse = " ")
    )
    })
    for (actCluster in 1:length(clusterFeatures)) {
    biclusteringSummary <- rbind(biclusteringSummary, clusterFeatures[[actCluster]])
    }
    colnames(biclusteringSummary) <-
    c(
      "DNDS columns",
      "Network columns",
      "Amount of Genes",
      "RAW DNDS",
      "STD DNDS",
      "0.1 Quantile DNDS",
      "0.9 Quantile DNDS",
      "Mean RAW DNDS",
      "STD RAW DNDS",
      "0.1 Quantile RAW DNDS",
      "0.9 RAW Quantile DNDS"
    )
    rownames(biclusteringSummary) <-
    paste0("Bicluster_", c(1:nrow(biclusteringSummary)))
    
    print(kable(biclusteringSummary,caption=gsub("_", " ", paste0(actBiclustering,"_",lowOrHigh))))

  }
}
ALL DNDS high
DNDS columns Network columns Amount of Genes RAW DNDS STD DNDS 0.1 Quantile DNDS 0.9 Quantile DNDS Mean RAW DNDS STD RAW DNDS 0.1 Quantile RAW DNDS 0.9 RAW Quantile DNDS
Bicluster_1 HUMAN LA-AVG 21 0.96 0.03 0.916 0.995 1.143 0.432 0.77 1.719
Bicluster_2 HUMAN WM-REST 19 0.962 0.03 0.915 0.995 1.178 0.44 0.768 1.745
Bicluster_3 HUMAN MO-RH+LH 18 0.959 0.026 0.926 0.989 1.065 0.377 0.789 1.478
Bicluster_4 HUMAN WM-REST WM-FACE 14 0.963 0.033 0.914 0.996 1.235 0.483 0.765 1.809
Bicluster_5 HUMAN cortico-limbic_network 13 0.954 0.028 0.926 0.992 1.04 0.366 0.787 1.62
Bicluster_6 MOCMNGPH_ND deafult_mode_network 8 0.949 0.036 0.912 0.988 0.945 0.196 0.758 1.176
Bicluster_7 MOUSE OTOGA salience_network 93 0.957 0.955 0.027 0.028 0.92 0.915 0.993 0.995 0.42 0.559 0.154 0.2 0.302 0.388 0.57 0.881
Bicluster_8 NEANDERTHAL GA-AVG 7 0.93 0.023 0.903 0.954 0.815 0.085 0.719 0.915
Bicluster_9 HUMAN GA-REWARD WM-REST 11 0.961 0.031 0.916 0.995 1.161 0.454 0.77 1.719
Bicluster_10 HUMAN EM-FACES MO-RH+LH 9 0.963 0.031 0.915 0.988 1.142 0.337 0.768 1.452
Bicluster_11 MOCMNGPH_ND sensorimotor_network 5 0.945 0.039 0.91 0.986 0.923 0.209 0.751 1.154
Bicluster_12 MOCMNDHN_PG visuospatial_network 4 0.958 0.03 0.935 0.987 1.274 0.556 0.915 1.797
Bicluster_13 OTOGA mPFC-AMY_network mPFC-Acb_network 31 0.944 0.031 0.91 0.989 0.507 0.174 0.377 0.718
Bicluster_14 MOCMNGPH_ND EM-FACES MO-RH+LH RE-AVG 4 0.948 0.035 0.923 0.982 1.023 0.419 0.796 1.408
Bicluster_15 DENISOVAN fronto-parietal 44 0.957 0.03 0.918 0.994 1.884 1.042 1.172 2.937
Bicluster_16 MOCMNGP_HND central_executive_network 27 0.948 0.034 0.907 0.996 1.393 0.711 0.834 2.549
Bicluster_17 DENISOVAN ventral_attention 37 0.954 0.026 0.923 0.99 1.663 0.568 1.198 2.479
Bicluster_18 NOMLE MO-RF+LF 26 0.947 0.029 0.906 0.98 0.869 0.266 0.597 1.15
Bicluster_19 DENISOVAN dorsal_attention 16 0.95 0.027 0.922 0.983 1.594 0.543 1.194 2.099
Bicluster_20 PANTR SO-TOM 7 0.958 0.036 0.914 0.993 1.446 0.7 0.85 2.154
ALL DNDS low
DNDS columns Network columns Amount of Genes RAW DNDS STD DNDS 0.1 Quantile DNDS 0.9 Quantile DNDS Mean RAW DNDS STD RAW DNDS 0.1 Quantile RAW DNDS 0.9 RAW Quantile DNDS
Bicluster_1 MOCMNGPD_HN LA-AVG 281 0.001 0.006 0 0 0.002 0.013 0 0
Bicluster_2 MOCMNGPD_HN WM-REST 222 0.001 0.008 0 0 0.002 0.016 0 0
Bicluster_3 MOCMNGPD_HN cortico-limbic_network 215 0.001 0.007 0 0 0.002 0.015 0 0
Bicluster_4 MOCMNGPD_HN deafult_mode_network 199 0 0.005 0 0 0.001 0.01 0 0
Bicluster_5 MOCMNGPD_HN salience_network 179 0.003 0.014 0 0 0.005 0.025 0 0
Bicluster_6 MOCMNGPH_ND WM-FACE 179 0.001 0.008 0 0 0.002 0.017 0 0
Bicluster_7 MOCMNGPH_ND GA-AVG 181 0.001 0.008 0 0 0.003 0.018 0 0
Bicluster_8 MOCMNGPH_ND MO-RH+LH 160 0.001 0.01 0 0 0.003 0.018 0 0
Bicluster_9 MOCMNGPH_ND EM-FACES 105 0.001 0.01 0 0 0.002 0.015 0 0
Bicluster_10 MOCMNGPD_HN GA-REWARD 73 0.002 0.012 0 0 0.005 0.024 0 0
Bicluster_11 MOCMNGPD_HN mPFC-AMY_network 66 0.002 0.011 0 0 0.004 0.021 0 0
Bicluster_12 MOCMNGPD_HN visuospatial_network 64 0 0 0 0 0 0 0 0
Bicluster_13 MOCMNGPD_HN sensorimotor_network 64 0 0 0 0 0 0 0 0
Bicluster_14 MOCMNGPD_HN fronto-parietal 62 0 0.004 0 0 0.002 0.013 0 0
Bicluster_15 MOCMNGPH_ND central_executive_network 70 0.002 0.011 0 0 0.005 0.023 0 0
Bicluster_16 MOCMNGPD_HN mPFC-Acb_network 54 0.001 0.008 0 0 0.002 0.016 0 0
Bicluster_17 MOCMNGPH_ND RE-AVG 51 0.002 0.013 0 0 0.003 0.02 0 0
Bicluster_18 MOCMNGPD_HN ventral_attention 37 0 0 0 0 0 0 0 0
Bicluster_19 NEANDERTHAL MO-RF+LF 31 0 0 0 0 0 0 0 0
Bicluster_20 MOCMNGPD_HN dorsal_attention 19 0 0 0 0 0 0 0 0
Bicluster_21 MOUSE OTOGA GORGO MO_CMNGPDNH MOC_MNGPDNH MOCMN_GPDNH LA-AVG MO-RH+LH MO-T 4 0.03 0.009 0 0 0 0 0.034 0.017 0 0 0 0.001 0.01 0 0 0 0 0 0.061 0.024 0 0 0 0.001 0.009 0.003 0 0 0 0 0.007 0.007 0 0 0 0 0.005 0 0 0 0 0 0.016 0.01 0 0 0 0
C mouse to human high
DNDS columns Network columns Amount of Genes RAW DNDS STD DNDS 0.1 Quantile DNDS 0.9 Quantile DNDS Mean RAW DNDS STD RAW DNDS 0.1 Quantile RAW DNDS 0.9 RAW Quantile DNDS
Bicluster_1 HUMAN LA-AVG 21 0.96 0.03 0.916 0.995 1.143 0.432 0.77 1.719
Bicluster_2 HUMAN WM-REST 19 0.962 0.03 0.915 0.995 1.178 0.44 0.768 1.745
Bicluster_3 HUMAN MO-RH+LH 18 0.959 0.026 0.926 0.989 1.065 0.377 0.789 1.478
Bicluster_4 HUMAN WM-REST WM-FACE 14 0.963 0.033 0.914 0.996 1.235 0.483 0.765 1.809
Bicluster_5 MOCMNGPH_ND cortico-limbic_network 8 0.955 0.033 0.912 0.988 0.967 0.182 0.758 1.176
Bicluster_6 MOCMNGPH_ND deafult_mode_network 8 0.949 0.036 0.912 0.988 0.945 0.196 0.758 1.176
Bicluster_7 MOCMN_GPDNH salience_network 129 0.949 0.028 0.91 0.985 1.25 1.172 0.713 1.653
Bicluster_8 HUMAN GA-AVG LA-AVG 11 0.959 0.026 0.94 0.992 1.065 0.356 0.817 1.619
Bicluster_9 HUMAN GA-REWARD WM-REST 11 0.961 0.031 0.916 0.995 1.161 0.454 0.77 1.719
Bicluster_10 HUMAN EM-FACES MO-RH+LH 9 0.963 0.031 0.915 0.988 1.142 0.337 0.768 1.452
Bicluster_11 MOCMNGPH_ND sensorimotor_network 5 0.945 0.039 0.91 0.986 0.923 0.209 0.751 1.154
Bicluster_12 MOCMNDHN_PG visuospatial_network 4 0.958 0.03 0.935 0.987 1.274 0.556 0.915 1.797
Bicluster_13 MOCMNGPH_ND EM-FACES MO-RH+LH RE-AVG 4 0.948 0.035 0.923 0.982 1.023 0.419 0.796 1.408
Bicluster_14 MOCMNGP_HND central_executive_network 27 0.948 0.034 0.907 0.996 1.393 0.711 0.834 2.549
Bicluster_15 MOC_MNGPDNH mPFC-Acb_network 37 0.943 0.03 0.905 0.983 1.109 0.692 0.668 1.622
Bicluster_16 MOC_MNGPDNH mPFC-AMY_network 33 0.948 0.032 0.91 0.993 1.284 0.902 0.691 2.356
Bicluster_17 MOC_MNGPDNH fronto-parietal 30 0.948 0.029 0.906 0.987 1.154 0.698 0.673 1.799
Bicluster_18 MOCMNGPH_ND LA-AVG cortico-limbic_network ventral_attention 3 0.954 0.019 0.938 0.966 0.915 0.094 0.837 0.983
Bicluster_19 MO_CMNGPDNH MO-RF+LF 24 0.953 0.029 0.912 0.987 0.694 0.308 0.451 0.895
Bicluster_20 MOCM_NGPDNH dorsal_attention 9 0.967 0.021 0.946 0.99 1.503 0.529 1.03 2.34
Bicluster_21 MOCM_NGPDNH MO-RH+LH SO-TOM WM-REST 5 0.952 0.04 0.912 0.992 2.698 3.594 0.78 6.144
C mouse to human low
DNDS columns Network columns Amount of Genes RAW DNDS STD DNDS 0.1 Quantile DNDS 0.9 Quantile DNDS Mean RAW DNDS STD RAW DNDS 0.1 Quantile RAW DNDS 0.9 RAW Quantile DNDS
Bicluster_1 MOCMNGPD_HN LA-AVG 281 0.001 0.006 0 0 0.002 0.013 0 0
Bicluster_2 MOCMNGPD_HN WM-REST 222 0.001 0.008 0 0 0.002 0.016 0 0
Bicluster_3 MOCMNGPD_HN cortico-limbic_network 215 0.001 0.007 0 0 0.002 0.015 0 0
Bicluster_4 MOCMNGPD_HN deafult_mode_network 199 0 0.005 0 0 0.001 0.01 0 0
Bicluster_5 MOCMNGPD_HN salience_network 179 0.003 0.014 0 0 0.005 0.025 0 0
Bicluster_6 MOCMNGPH_ND WM-FACE 179 0.001 0.008 0 0 0.002 0.017 0 0
Bicluster_7 MOCMNGPH_ND GA-AVG 181 0.001 0.008 0 0 0.003 0.018 0 0
Bicluster_8 MOCMNGPH_ND MO-RH+LH 160 0.001 0.01 0 0 0.003 0.018 0 0
Bicluster_9 MOCMNGPH_ND EM-FACES 105 0.001 0.01 0 0 0.002 0.015 0 0
Bicluster_10 MOCMNGPD_HN GA-REWARD 73 0.002 0.012 0 0 0.005 0.024 0 0
Bicluster_11 MOCMNGPD_HN mPFC-AMY_network 66 0.002 0.011 0 0 0.004 0.021 0 0
Bicluster_12 MOCMNGPH_ND visuospatial_network 86 0.002 0.01 0 0 0.004 0.021 0 0
Bicluster_13 MOCMNGPD_HN sensorimotor_network 64 0 0 0 0 0 0 0 0
Bicluster_14 MOCMNGPD_HN fronto-parietal 62 0 0.004 0 0 0.002 0.013 0 0
Bicluster_15 MOCMNGPD_HN central_executive_network 58 0.003 0.013 0 0 0.006 0.026 0 0
Bicluster_16 MOCMNGPD_HN mPFC-Acb_network 54 0.001 0.008 0 0 0.002 0.016 0 0
Bicluster_17 MOCMNGPH_ND RE-AVG 51 0.002 0.013 0 0 0.003 0.02 0 0
Bicluster_18 MOCMNGPD_HN ventral_attention 37 0 0 0 0 0 0 0 0
Bicluster_19 MOCMNGPH_ND MO-RF+LF 34 0 0 0 0 0 0 0 0
Bicluster_20 MOCMNGPD_HN dorsal_attention 19 0 0 0 0 0 0 0 0
Bicluster_21 MO_CMNGPDNH MOC_MNGPDNH MOCM_NGPDNH MOCMN_GPDNH LA-AVG MO-RH+LH MO-T 4 0 0 0 0 0 0 0 0.001 0 0 0 0 0 0 0 0.001 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
E2 hominids speciation with ancestorsNOG high
DNDS columns Network columns Amount of Genes RAW DNDS STD DNDS 0.1 Quantile DNDS 0.9 Quantile DNDS Mean RAW DNDS STD RAW DNDS 0.1 Quantile RAW DNDS 0.9 RAW Quantile DNDS
Bicluster_1 HUMAN LA-AVG 21 0.96 0.03 0.916 0.995 1.143 0.432 0.77 1.719
Bicluster_2 NEANDERTHAL WM-REST 10 0.939 0.026 0.907 0.962 0.87 0.145 0.735 0.957
Bicluster_3 HUMAN MO-RH+LH 18 0.959 0.026 0.926 0.989 1.065 0.377 0.789 1.478
Bicluster_4 DENISOVAN salience_network 188 0.946 0.03 0.907 0.99 1.634 0.742 1.122 2.518
Bicluster_5 DENISOVAN cortico-limbic_network deafult_mode_network 106 0.945 0.027 0.91 0.983 1.543 0.584 1.14 2.07
Bicluster_6 NEANDERTHAL GA-AVG 7 0.93 0.023 0.903 0.954 0.815 0.085 0.719 0.915
Bicluster_7 HUMAN GA-REWARD WM-REST 11 0.961 0.031 0.916 0.995 1.161 0.454 0.77 1.719
Bicluster_8 HUMAN EM-FACES MO-RH+LH 9 0.963 0.031 0.915 0.988 1.142 0.337 0.768 1.452
Bicluster_9 DENISOVAN sensorimotor_network 62 0.951 0.028 0.916 0.989 1.731 0.901 1.164 2.442
Bicluster_10 DENISOVAN mPFC-Acb_network 45 0.951 0.026 0.919 0.984 1.628 0.742 1.176 2.123
Bicluster_11 MOCMNGPH_ND EM-FACES MO-RH+LH RE-AVG 4 0.948 0.035 0.923 0.982 1.023 0.419 0.796 1.408
Bicluster_12 DENISOVAN fronto-parietal 44 0.957 0.03 0.918 0.994 1.884 1.042 1.172 2.937
Bicluster_13 MOCMNGP_HND central_executive_network 27 0.948 0.034 0.907 0.996 1.393 0.711 0.834 2.549
Bicluster_14 DENISOVAN mPFC-AMY_network 39 0.944 0.028 0.908 0.981 1.561 0.788 1.126 2.029
Bicluster_15 DENISOVAN ventral_attention 37 0.954 0.026 0.923 0.99 1.663 0.568 1.198 2.479
Bicluster_16 DENISOVAN MO-RF+LF MO-RH+LH 27 0.953 0.028 0.915 0.992 1.672 0.586 1.157 2.641
Bicluster_17 DENISOVAN dorsal_attention 16 0.95 0.027 0.922 0.983 1.594 0.543 1.194 2.099
Bicluster_18 PANTR SO-TOM 7 0.958 0.036 0.914 0.993 1.446 0.7 0.85 2.154
E2 hominids speciation with ancestorsNOG low
DNDS columns Network columns Amount of Genes RAW DNDS STD DNDS 0.1 Quantile DNDS 0.9 Quantile DNDS Mean RAW DNDS STD RAW DNDS 0.1 Quantile RAW DNDS 0.9 RAW Quantile DNDS
Bicluster_1 MOCMNGPD_HN LA-AVG 281 0.001 0.006 0 0 0.002 0.013 0 0
Bicluster_2 MOCMNGPD_HN WM-REST 222 0.001 0.008 0 0 0.002 0.016 0 0
Bicluster_3 MOCMNGPD_HN cortico-limbic_network 215 0.001 0.007 0 0 0.002 0.015 0 0
Bicluster_4 NEANDERTHAL deafult_mode_network 190 0.004 0.017 0 0 0.005 0.025 0 0
Bicluster_5 MOCMNGPD_HN salience_network 179 0.003 0.014 0 0 0.005 0.025 0 0
Bicluster_6 MOCMNGPD_HN WM-FACE 134 0.001 0.006 0 0 0.002 0.014 0 0
Bicluster_7 NEANDERTHAL GA-AVG 122 0.003 0.014 0 0 0.006 0.024 0 0
Bicluster_8 NEANDERTHAL MO-RH+LH 115 0.002 0.011 0 0 0.003 0.018 0 0
Bicluster_9 NEANDERTHAL EM-FACES 81 0 0 0 0 0 0 0 0
Bicluster_10 MOCMNGPD_HN GA-REWARD 73 0.002 0.012 0 0 0.005 0.024 0 0
Bicluster_11 MOCMNGPD_HN mPFC-AMY_network 66 0.002 0.011 0 0 0.004 0.021 0 0
Bicluster_12 MOCMNGPD_HN visuospatial_network 64 0 0 0 0 0 0 0 0
Bicluster_13 MOCMNGPD_HN sensorimotor_network 64 0 0 0 0 0 0 0 0
Bicluster_14 NEANDERTHAL central_executive_network 56 0.001 0.007 0 0 0.003 0.015 0 0
Bicluster_15 MOCMNGPD_HN fronto-parietal 62 0 0.004 0 0 0.002 0.013 0 0
Bicluster_16 MOCMNGPD_HN mPFC-Acb_network 54 0.001 0.008 0 0 0.002 0.016 0 0
Bicluster_17 MOCMNGPD_HN RE-AVG 40 0 0 0 0 0 0 0 0
Bicluster_18 MOCMNGPD_HN ventral_attention 37 0 0 0 0 0 0 0 0
Bicluster_19 NEANDERTHAL MO-RF+LF 31 0 0 0 0 0 0 0 0
Bicluster_20 MOCMNGPD_HN dorsal_attention 19 0 0 0 0 0 0 0 0
C3 mouse to human late no gorgo high
DNDS columns Network columns Amount of Genes RAW DNDS STD DNDS 0.1 Quantile DNDS 0.9 Quantile DNDS Mean RAW DNDS STD RAW DNDS 0.1 Quantile RAW DNDS 0.9 RAW Quantile DNDS
Bicluster_1 MOCMNGP_HND LA-AVG 102 0.944 0.029 0.908 0.987 1.247 0.551 0.838 1.893
Bicluster_2 MOCMNGP_HND cortico-limbic_network 85 0.948 0.031 0.909 0.995 1.354 0.647 0.846 2.353
Bicluster_3 MOCMNGP_HND WM-REST 79 0.947 0.029 0.908 0.985 1.283 0.554 0.837 1.767
Bicluster_4 HUMAN MO-RH+LH 18 0.959 0.026 0.926 0.989 1.065 0.377 0.789 1.478
Bicluster_5 MOCMNGP_HND deafult_mode_network 72 0.948 0.031 0.91 0.992 1.327 0.612 0.854 2.117
Bicluster_6 MOCMNGP_HND salience_network 71 0.953 0.026 0.921 0.989 1.313 0.44 0.915 1.982
Bicluster_7 MOCMNGP_HND WM-FACE 54 0.951 0.027 0.914 0.984 1.326 0.579 0.876 1.715
Bicluster_8 HUMAN GA-REWARD WM-REST 11 0.961 0.031 0.916 0.995 1.161 0.454 0.77 1.719
Bicluster_9 HUMAN GA-AVG LA-AVG 11 0.959 0.026 0.94 0.992 1.065 0.356 0.817 1.619
Bicluster_10 MOCMNGP_HND EM-FACES 34 0.956 0.027 0.917 0.984 1.349 0.481 0.896 1.728
Bicluster_11 MOCMNGP_HND central_executive_network 27 0.948 0.034 0.907 0.996 1.393 0.711 0.834 2.549
Bicluster_12 MOCMNGP_HND visuospatial_network 27 0.948 0.032 0.904 0.986 1.265 0.414 0.821 1.833
Bicluster_13 MOCMNGP_HND sensorimotor_network 24 0.943 0.029 0.904 0.985 1.215 0.51 0.824 1.832
Bicluster_14 MOCMNGP_HND mPFC-AMY_network 22 0.945 0.028 0.913 0.988 1.297 0.71 0.874 1.92
Bicluster_15 MOCMNGP_HND mPFC-Acb_network 20 0.952 0.028 0.916 0.995 1.407 0.743 0.887 2.42
Bicluster_16 MOCMNGP_HND fronto-parietal 19 0.948 0.03 0.907 0.981 1.295 0.546 0.837 1.736
Bicluster_17 MOCMNGP_HND RE-AVG 19 0.956 0.033 0.913 0.994 1.484 0.685 0.868 2.299
Bicluster_18 MOCMNGP_HND MO-RF+LF 16 0.962 0.029 0.926 0.994 1.578 0.779 0.942 2.662
Bicluster_19 MOCMNGPH_ND LA-AVG cortico-limbic_network ventral_attention 3 0.954 0.019 0.938 0.966 0.915 0.094 0.837 0.983
Bicluster_20 MOCMNGP_HND SO-TOM WM-REST 3 0.933 0.03 0.911 0.958 1.033 0.256 0.858 1.245
Bicluster_21 MOCMNGP_HND fronto-parietal dorsal_attention 3 0.953 0.031 0.927 0.973 1.248 0.318 0.983 1.478
C3 mouse to human late no gorgo low
DNDS columns Network columns Amount of Genes RAW DNDS STD DNDS 0.1 Quantile DNDS 0.9 Quantile DNDS Mean RAW DNDS STD RAW DNDS 0.1 Quantile RAW DNDS 0.9 RAW Quantile DNDS
Bicluster_1 HUMAN LA-AVG 493 0.001 0.008 0 0 0.003 0.016 0 0
Bicluster_2 HUMAN WM-REST 395 0.001 0.006 0 0 0.002 0.014 0 0
Bicluster_3 HUMAN cortico-limbic_network 357 0.002 0.01 0 0 0.003 0.017 0 0
Bicluster_4 HUMAN deafult_mode_network 341 0.001 0.01 0 0 0.003 0.017 0 0
Bicluster_5 HUMAN salience_network 292 0.002 0.01 0 0 0.003 0.018 0 0
Bicluster_6 HUMAN WM-FACE 248 0 0.004 0 0 0.001 0.011 0 0
Bicluster_7 HUMAN GA-AVG 242 0.001 0.007 0 0 0.003 0.017 0 0
Bicluster_8 HUMAN MO-RH+LH 229 0.001 0.005 0 0 0.002 0.012 0 0
Bicluster_9 HUMAN EM-FACES 144 0.001 0.005 0 0 0.002 0.013 0 0
Bicluster_10 HUMAN GA-REWARD 129 0.002 0.008 0 0 0.004 0.019 0 0
Bicluster_11 HUMAN visuospatial_network 114 0.002 0.012 0 0 0.005 0.022 0 0
Bicluster_12 HUMAN sensorimotor_network 110 0.002 0.012 0 0 0.005 0.022 0 0
Bicluster_13 HUMAN central_executive_network 99 0.001 0.007 0 0 0.001 0.012 0 0
Bicluster_14 HUMAN mPFC-AMY_network 91 0.002 0.014 0 0 0.003 0.019 0 0
Bicluster_15 HUMAN mPFC-Acb_network 82 0.002 0.014 0 0 0.003 0.021 0 0
Bicluster_16 HUMAN RE-AVG 73 0.001 0.004 0 0 0.002 0.014 0 0
Bicluster_17 HUMAN ventral_attention 60 0.001 0.009 0 0 0.002 0.015 0 0
Bicluster_18 HUMAN MO-RF+LF 46 0 0 0 0 0 0 0 0
Bicluster_19 HUMAN fronto-parietal dorsal_attention 28 0 0 0 0 0 0 0 0
Bicluster_20 MOCMNGP_HND SO-TOM 37 0.007 0.017 0 0.029 0.009 0.021 0 0.048